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Notation 


No NU {0}. 

f The Fourier transform of the function f. 

fe The adjoint/complex conjugate of the function/value f. 
MB The Hermitian conjugate of the vector/matrix M. 

F The discrete Fourier transform of the discrete function F. 
(A, B] The Lie bracket or commutator |A, B] = AB — BA. 


[A,B,C,D] The repeated Lie bracket or commutator with the form [A, [B, [C, D]]]. 


1: Introduction and Preliminary 


Mathematics 


1.1 Introduction 


Splitting methods for ordinary and partial differential equations (PDEs and ODEs) have 
a long history although, it is only recently that they have appeared under this name and 
that the theory has been assembled into the unified field that exists today. Splitting meth- 
ods were originally developed for the traditional numerical motivations of speed, accuracy, 
and stability. It is now clear, however, that splitting methods are also an excellent general 
and flexible way of constructing geometric integrators; integrators for vectors fields which 


posses certain geometric properties, such as being Hamiltonian or conserving certain quan- 


tities, which one wishes to preserve |McLachlan and Quispel, 2001 [Budd and Iserles, 1999 
Hairer et al., 2002|. By using splitting methods to construct integrators which preserve 


structural features of a flow X one is able to confer superior qualitative properties to an 
integrator, especially when integrating for long times. This has led to splitting methods 
becoming the method of choice in a variety of fields from celestial mechanics to molecular 
dynamics and particle accelerator physics. Splitting methods can yield integrators which, 


on grounds of their efficient implementation, are attractive even when one ignores their 


geometric properties. McLachlan |McLachlan and Quispel, 2001| gives the following caveat: 


“Because we’re asking something more of our method, it may turn out to be 


computationally more expensive than a standard method.” 
but then follows this warning with the providential statement 


“Amazingly, (because the laws are so natural?), sometimes it’s actually much 


cheaper.” 
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1.1.1 What is splitting? 


Splitting methods are concerned with the time integration of ODEs and PDEs. For a phase 
space M, differential equation t = X(x), r € M, and X a vector field on M, splitting 


methods involve three equally important steps: 
1. choosing a set of vectors fields X; such that we can decompose X as X = У? Xj; 
2. integrating, either exactly or approximately, each of the X;; and 
3. combining these solutions to yield an integrator for X. 


An example: by writing the flow (i.e. the exact solution) of the ODE тї = X as z(t) = 


exp(t.X)z(0) we may use the composition method 
p(T) = ехр(тХу) exp(r X2) --- exp(r X4), 
where 7 is the time step, to approximate the flow. This is first order accurate in that 


(r) = exp(r У? X;) + O(r?). 


When splitting the vector field X, the pieces X; should always be chosen such that they 


are simpler than the original vector field. This can occur in two possible senses: 


1. The Xj; are of a simpler type than X. For example the Navier-Stokes equation 


д 1 
E = -uVu — A + Vu 


contains advection (Vu), pressure (V P), and diffusion (V?u) terms, each with distinct 


characteristic properties and appropriate numerical methods. 


2. The X; are the same type as X, but are easier to treat numerically. Examples are 
dimensional splitting for the multidimensional heat equation, Hamiltonian splitting for 
Hamiltonian ODEs and, importantly for us, the (linear) Schródinger equation iu, — 
Ure + V (x)u — each piece of which is linear and Hamiltonian, but the first term can be 


integrated more quickly, (in Fourier space, using the FFT), than the combined system. 


1.1.2 Historical development of splitting methods 


Although it is not a new field, the historical progress of splitting methods is far from clear 


with different developments occurring independently in various areas. Splitting methods 
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have developed in response to a variety of mathematical and physically motived stimuli 
with dimensional splitting for parabolic PDEs such as the multidimensional heat equa- 
tion, fractional-step and operator splitting methods for the Navier-Stokes equations and 
for reaction-diffusion systems, split-step methods in optics and acoustics, split Hamiltonian 
methods in chemical physics, the mapping method in celestial mechanics and Lie-Trotter- 
Kato formulae in quantum statistical mechanics (McLachlan and Quispel, 2002]. The field of 
splitting methods as it exists today is comprised of a combination of, and cross-over between 


all these areas. 


Splitting methods essentially began with the product formula of Trotter [Irotter, 1959 


lim (Ree — e AB). 
n=% 
where A and B are linear operators on a Banach space and either t € R or t € iR, t > 0, and A 
and B are bounded above. For parabolic equations key early references come from Tappert 
who introduced the split-step method for the 
nonlinear Schrödinger equation, the first example of what is now known as a symplectic 
integrator for a PDE. Another important idea in the development of splitting methods is the 
leapfrog or Störmer-Verlet method for Hamiltonian ODEs of the form 2 = f(x) = —VV (x). 
Although this method is usually attributed to Verlet who used it to simulate 
the molecular dynamics of Argon atoms interacting via a Lennard-Jones potential, (V (x) = 
is glo/rg) — (c/rijó, rij = ||; — |), the method can also be traced to Delambre 
and Newton [Newton, 1687]. A key development in splitting methods from 


within the physics community was the study of so-called “standard-mappings” which inspired 
the creation of the first symplectic integrators for celestial mechanics [Wisdom, 1982]. 


An important step in the development of splitting methods was the derivation of the 
leapfrog method, (already known to be a symplectic integrator and in common use as 
such at this point), as a composition of flows of elementary Hamiltonians, that is, as a 
splitting method. The systematic development of splitting methods by numerical ana- 
lysts was triggered by a preprint from Neri who used the Baker-Campbell- 
Hausdorff formula to derive composition methods of higher order. This was the inspira- 
tion for a more systematic development by Yoshida [Yoshida, 1990]. As an indication of 
both the growth in the field of splitting methods and of the significance of the publica- 
tion we point out that 15 years later has received well over 400 citations 
across diverse fields. In a parallel development, Suzuki began studying vari- 


ations and generalisations of the Trotter formula which led to new composition methods 


suzuki, 1990| Suzuki, 1992| Suzuki, 1994| Suzuki, 1995]. 
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1.1.3 Background to nonlinear oscillatory problems 


Nonlinear oscillatory problems often arise in the study of the propagation of waves, for 
example, gravity waves travelling in deep water or optical waves travelling in a dielectric. 
Such problems are often described by a PDE which can be split into linear and nonlinear 
parts and are therefore amenable to solution via splitting methods which deal separately 
with each of these parts. Much of the work in this area has been motivated by the study of 
optical waves propagating through media with a nonlinear response to the amplitude of the 
wave. An important example is the modelling of light waves in optical fibres for the simula- 
tion of optical communications systems 
[Sinkin et al., 2003]. Although the study of wave propagation through nonlinear oscillatory 


problems can be viewed from a number of contexts, we will focus on a description arising 


from the propagation of light waves through nonlinear media. 


Maxwell’s equations and wave propagation 


A theory for the propagation of light as a unified treatment of electric and magnetic fields 
was first presented by James Clerk Maxwell in the 1890s. His theory is given by a set of 
four equations relating the electric and magnetic intensity fields E and H, and the electric 


and magnetic induction fields D and B. In MKS шай Maxwell’s equation are 


М.Р = p, (1.1) 
Voices. ad (1.2) 
OB 
OD. 


where p and j are the free electric charge and current density respectively. In three- 


dimensional Cartesian coordinates the divergence and curl operators can be given by 








V.F= c 
Ox Oy Oz 
and 
д д 
— д д 
Vx Е = == 0 -— | F 
9 8 0 
Oy Ox 


'MKS is the system of units based on measuring length in metres, mass in kilogrammes and time in 
seconds. 
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for the vector field F = (Fy, Fy, Fz). 

For a dielectridd both p and j are zero. The four quantities E, H, D, and B are 
related by the so-called constitutive relations B = ИН and D = eg E + Р = cE. For 
optical waves in a dielectric it is normally assumed that the medium is non-magnetic and 
therefore u = uo, where ио is constant] [Newell and Moloney, 1992]. eis a dielectric constant 
which depends on the polarisation properties of the medium through which the light is 
propagating. More details on the constitutive relations can be found in chapter two of 
Butcher and Cotter, 1990}. 

When light propagates in a dielectric the electric field interacts with the medium through 
which it is travelling. This distorts the structure of the medium and induces a polarisation 
field P which depends on E. It is the nature of this dependence which is important in 
determining many of the properties concerning how the wave propagates. In particular, 
whether P depends on F in a linear or a nonlinear fashion. 

For small amplitudes of E the dependence is generally linear and takes the form P = egxE 
where x is known as the electric susceptibility coefficient. 

In order to derive a wave equation for the propagation of an electric field from Maxwell’s 
equations we take the curl of equation (L3) and change the order of the space and time 


derivatives on the right-hand term to get 


VERE 7 xB). 


ôt 


We then use the fact that В = oH to get 


Vx Vx B= -mŽ(V x H). 


We can now use the fact that j = 0 along with equation (LA) to replace the V x H term 
and get 
2 


VxVxE= D. (1.5) 


og 
The constitutive equation D = є0Е + P allows us to rewrite equation (L.L.3) as 


ФЕ PP 
VxVx E + cotto — —Н0- 

2А substance in which the valence and conduction bands are well separated and hence which is neither a 
conductor, nor a semi-conductor. 

З ро is known as the permeability of free space. In MKS units it has the value uo = 4r x 1077 = 
1.2566 x 10-°NA~? where N stands for Newtons and A for Amps. ио is related to the constant eo, known 
as the permittivity of free space, through the relationship eo — Tu where c = 3x 105m/ s is the speed of 
light in a vacuum. 
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which using the relationship jug = 1/c?eo becomes 


1 E 1 &P 
Е + = с. 1. 
VxVX Ets = As 58 (1.6) 


The left-hand-side of equation { Т.З) can be rewritten using the vector calculus identity 
VxVx E = V(V- E) — V?E (where V is the gradient operator and hence V? is the 


Laplacian which is sometimes written as A) to give 


180?E 1 Ө?Р 
v? V(V = 
= e o Dodge c?eg Ot? ` en 


This is the central equation in the description of electric field propagation in both linear and 


nonlinear media. 


For (isotropic, homogeneous) linear media the polarisation field P is given by 
P = єхЕ (1.8) 


where y is the linear electric susceptibility coefficient which may depend on frequency but 


which is independent of the intensity of the electric field. Substituting equation (L8) into 


(L7) gives 


Since x is independent of E equation (LI) of Maxwell’s’ equations along with p = 0 gives 
V-D=V. (ëE + P)= еу: (E - XE) = 0, 


and therefore V - E = —xV - E. However, on physical нш x > 0 and hence, for linear 
media V - E = 0. Equation (L7) then simplifies to 
n? OE 

VE-——— =0 1.9 

c2 Ot? ? ( ) 

where n = y1+x > 1 is the refractive index of the medium. In this case the electric field 


propagates through the medium with a phase velocity of c/n. 


‘If we allow x < 0 then the refractive index of a medium given by n = yI F x can be less than one. This 
implies that in such a medium with x < 0 the electric field travels faster than the speed of light in a vacuum. 
This would contradict Einstein’s famous theory of special relativity. 
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The dispersion relation and frequency dependence of the refractive index 


Since matter does not respond instantly to stimulation from the electric field the description 
of the polarisation should include some sort of delay term. For an isotropic medium this 
effect can be included through the equation Р = eo JTS x(t — s)E(s)ds where, due to the 
causality principle, x(t — s) = 0 for t — s < 0. The time dependence of the polarisation field 
and the susceptibility coefficient implies that that the Fourier transform of x(t) is frequency 
dependent (and complex). Details on the origins and the implications of this are discussed 
in (pp. 26-29). The main result of interest here is that the real 
part of X(w) means that different frequencies of light propagate at different speeds. 

The most basic solution, from which other solutions of the wave equation can be con- 


structed, is a singly polarised plane wave or wavetrain. In one dimension this takes the form 


E(a,t) = Aexp(i(ka — wt)) + c.c., (1.10) 


where A € C, k is the wave number (or wave vector in the case of more than one space 
dimension) and c.c. denotes the complex conjugate of the preceding term. Substituting 
(10) into equation (LI) it is easy to see that in order to be a solution { ТО) must satisfy 


d(k,w)A = (e — ише) А = 0. (1.11) 


Since A Z 0 we must therefore have d(k,w) = 0. In general n(w) is complex but in most cases 
the imaginary component is many orders of magnitude smaller than the real component and 
can be ignored. We therefore assume here that n(w) is a real valued function. 

The relationship d(k,w) = 0 is known as a dispersion relation. For nonlinear media the 
dispersion relation also depends on the square of the amplitude A of the electric field. For 
plane wave solutions d = 0 is algebraic. We will deal with this case in chapter [2] 

When the solution of (L7) is not a plane wave, but rather a wave packet, the relationship 
d — 0 takes the form of a PDE which describes the envelope of the wave packet. When the 
medium in question is nonlinear, then so too is the PDE. We will treat this case in more 
detail in chapter [5] 


'The relationship between P and E for nonlinear media 


When the electric field in a medium is intense, the electric susceptibility x, and hence also 
the refractive index n of the medium depend not only on the frequency, but also on the 
intensity I of the electric field. Since J = |A|?, where A is the amplitude of the electric field, 


this gives a nonlinear relationship between n and A. As a result, the intensity of the wave 
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has a feed-back effect and can alter the phase velocity of the wave. Physical causes which 
give rise to the nonlinear refractive index include the electronic and vibrational resonances 
of atoms and molecules as well as changes in density due to heating of the medium. 

For nonlinear media a relationship between P and E can be found by writing P as the 
sum of a linear Pr, and a nonlinear Pyr, part. In the one-dimensional case the equation for 


propagation of waves in a nonlinear Kerr sedit equation (LJ) becomes 


O^E(rt)) 16°E(a2,t) 1899 Г ax 1 а? 
Ox? 2 Ot? c2 Ot? 


t—s)E(a,s)ds = adag ҮР (1.12) 
where Pr is given by the integral term and where y“!) is the first order susceptibility tensor. 
To leading order the Pyz term depends on the higher order susceptibility tensor x®). The 
properties of these tensors are determined by the physical structure of the medium through 
which the wave is propagating, in particular, by any symmetry within the structure. Details 
concerning these properties can be found in chapter two of (Shen, 1984], chapter four of 
and chapters one and three of [Boyd, 1992]. For a Kerr medium 


one assumes that |x) [EP « x® |El. 


It is shown in [Newell and Moloney, 1992], (pp.34-36), that from equation (L.I2) one 


arrives at a nonlinear expression for the refractive index of a medium. This is given by 
n(w, AP) = no(w) + na(w)|AP, — 2nona = 3° (w), (1.13) 


where no(w) and na(w) are the refractive index coefficients of order zero and two respectively, 
A is the amplitude of E and $9 (w) is the Fourier transform of x 9). Substituting equation 
(L.13) into the dispersion relation (LIT) gives the following nonlinear dispersion relation, 


[n Co, [А|?)]#ш° 


d(k,w,|A|?) = k? — - =0. (1.14) 


C 


Wave packets 


An important concept in the study of wave propagation in both linear and nonlinear media 
is the idea of wave packets. Light waves are rarely exactly monochromatic, (i.e. of only a 
single frequency), and more usually consist of a linear combination of plane waves whose 
wave numbers, (or wave vectors), and frequencies lie in a narrow band of order 0 < u & 1 


about some central values kg and wo. Such a collection is commonly called a wavepacket 


5A Kerr medium is one in which the nonlinear response of the medium is small with respect to the linear 
response. 
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and is represented in one-dimension by 
E(x,t) = A(x, t) ехр( (Кох — wot)) + с.с. + О(и), d(ko,wo) = 0, (1.15) 


where the complex envelope A(x, Ё) is a slowly varying function of space and time. That is, 
OA/Ox and OA/Ot are of order u|ko| A and и|шо|А respectively and hence, 


0?A 
OU 


0?A 
Oa? 








«x « |A], 








«m « |A]. 


To see this we represent the electric field as a linear combination of plane waves with wave 
numbers and frequencies ko + ШЕ; and wo + uQ; such that the dispersion relation given by 
(LIT) is still satisfied. This gives 


E(x,t) = > A; exp(i(ko + uK;)x — (wo + uN;)t) + с.с. 
j 
which can be rewritten as (1.15) by setting A(z,t) = У, Aj exp(i(uKjx — uNjt)). That is, 


A(z,t) is a function of uz and pt and, therefore, a slowly varying function of x and t. 


In chapter BJ] we will use the points mentioned here about wavepackets, dispersion and 
nonlinear wavetrains in Kerr media to derive an equation for the propagation of optical 
pulses in nonlinear media which includes second order dispersion effects. We then show how 
splitting methods can be used to find numerical solutions of this nonlinear wave equation. 
In chapter [2] we use a first order splitting method and the nonlinear Schródinger equation 
(NLSE) to illustrate how splitting methods can be applied to nonlinear problems where the 
PDE can be split into a linear and a nonlinear term. We consider two different methods for 
the spatial discretisation of the problem; a finite-difference method and a pseudo-spectral 
Fourier method. Both of these methods are analysed with respect to stability, conserva- 
tion and dispersion and the results compared with the analytic properties of the nonlinear 
Schrödinger equation. Numerical results are given for plane wave solutions. In chapter [8] we 
focus on the split-step Fourier method and extend our numerical results to include soliton 
solutions of the NLSE. Higher order splitting methods are developed in chapter Й and im- 
plemented as split-step Fourier methods. These methods are again applied to the nonlinear 
Schródinger equation and the numerical results are compared with an exact soliton solution 
in order to illustrate some of the overall behaviour of the methods. We use the remainder 


of this chapter to present some mathematical results which we will later find useful. 
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1.2 Preliminary mathematics 


1.2.1 The discrete Fourier transform 


Assume we have a function f which is defined on R and which is absolutely integrable, i.e. 
[I Mf (2)|dx < oo. We can then define a function 


Ўш) = ГИ f (x) exp(—i2mwz)dx (1.16) 


with —oo < ш < oo. The function (LI6) is the (unique) Fourier transform of f(x). We 
will refer to j (w) as being defined in the frequency domain and f(x) as being defined in 
the spatial domain, or in the temporal domain when f is a time dependent function. The 


function f can be recovered from f via the inverse Fourier transform, 
eo ^ 
f(x) x f (w) xpli2rwe)dw. (1.17) 
00 


If, additionally, f is a continuous, L—-periodic function, і.е. f(x + L) = f(x), n € Z then we 


can define f(x) as a Fourier series, 


оо 


fo) > сьехр(іџьх) (1.18) 
k=—0o 
where a 
= i. f(x) exp(—ipzx)dx (1.19) 


and ик = 276/1. If f(x) is only piecewise continuous, but the one-sided limits 
7)= li ith) and (а) = li ith 
Flay) = lim f(z-Fh) and Ла) = lim f(a +h) 


$ " E f xt +f xr 
exist for all x; then we must replace the left-hand side of (LIS) with foro) 


The interval [- L/2, L/2] can be discretised by dividing it into N equal sub-intervals 
B ana ж = nh for n = —N/2,...,N/2. 
After setting g(x) = f(x) exp(-iuxx) the integral in (I9) can be approximated by the 


with a grid spacing h = L/N, where N is eve 


6 Although we assume here that N is always even (and positive), all the results also hold for N odd if the 
expressions (N/2) — 1 are changed to (N — 1)/2. 
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trapezoidal rule. This gives 





1, А1 
1 fa 1А -L \ L 
Ck — T n К g(x)dz ~ 1519 (+) +2 У; g(tn) +g (5) (1.20) 
2 n=Z +1 
If f is only piecewise continuous with a possible discontinuity at 2 = +L/2, then when the 
Fourier series converges, the values at the endpoints 2 = +L/2 converge to the average value 








iC) 


Since exp(ipja) = (—1)* when x = +L/2 the function g(x) is altered to read 





б) = f (x) &p(-iurz) for x z +4 (L3) 
nz c» (f (=) +/ (&)) for x = +3 | 


which guarantees that g(-L/2) = g(L/2). The function g(x) must be treated similarly 





at any internal grid points where a discontinuity exists. Using equation (1.21) in equation 


(20) gives the following expression for the Fourier coefficients 


Ny Nj 
1, x lx 
Ck < т^ у е N у lan) xpl-iukEn). (1.22) 
—N -N 
сор NE 


The discrete Fourier transform and its inverse 


Comparing equations (L.16) and (£19) shows that if f(x) is a spatially limited function 
such that f(x) = 0 for |x| > L/2, then the Fourier transform of f(x), evaluated at the 
frequency wy = k/L, is given by for) = Гек. For the discrete function F defined by 
Fn = f(an), n = —N/2,..., №2 — 1 the discrete Fourier transform (DFT) is then defined 
by 

X 


`> Fn ехр(—ішьхъ), k = =g ete 
-N 


n= 


^ 1 


№ = = —1. (1.23) 


and hence f(wr) > LF. The inverse discrete Fourier transform (IDFT) is defined analo- 


gously as 


Fn = > Fy, ехр(—їшһть), M= p 


LIN 
2 


-1 (1.24) 


> 
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If the Fourier series given by (LIS) is evaluated at the points x, and is used in the 
calculation of the DFT (1.23) then 


2 oo 
^ 1 


Fy = E > > Cj exp(ipj2n) exp(—tpUp2n) 


LIN j=—-oo 
n= 2 


Il 
2 
=|е 
Ф 
Ы 
s 
1 
bo 
A 
zm 
| 
3 
3 


[ 
М 
= 
3 
= 


(1.25) 


where the last equality follows from the fact that, for (k — j) Z 0, exp(—i2a(k — j)n/N) are 
the N—th roots of unity and hence 


NO 
у; =. =й, N when kz j mod N 
exp ——————) = 
=" N 0 otherwise. 
n= 


From this result, one can show that if F is the sequence obtained by periodically sampling 
f(x), and if f(x) is L-periodic with p — 1 continuous derivatives and a piecewise continuous 
p—th derivative, then Ё, = f (an) + O(hP). 


Parseval’s equality 


Given two, (possibly complex valued), functions f(x) and g(x), and their Fourier transforms 


f(w) and g(w), Parseval’s equality states 


f "Pooks f "Ponai (1.26) 


where f* denotes the complex conjugate of the function f. The equality is easily verified; 


f ro Cf оерт) а 
= [ao (f FO erorar) az 
ІА 2x (. fe) exp(-i2mutr)dr) 7 


f аслас)" 


00 


ГА _ /"(ж)д(®)й 
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If we set g(x) = f(x) then equation (26) gives || f(x)||? = || f(w)||2 where ||: || denotes the 
L?-norm. If the continuous function f(x) is replaced with its discrete analogue F as defined 


earlier then the discrete version of Parseval’s equality is given by 


N- N 1 

2 1 2 76 

> 15е У Al. (1.27) 
а= а= 


1.2.2 Тһе Baker-Campbell-Hausdorff formula 


Given two arbitrary, generally non-commuting matrices А and В we want to find a matrix 
C(t) such that for t real and sufficiently small 


exp(tA) exp(tB) = exp(C(t)). (1.28) 


where the matrix exponential is given by the Taylor dere 


exp(A) = ur 
k=0 ` 


In order to get a first idea of the general form of C(t) we expand the left-hand side of (28) 
as a product of two such series: 

exp(tA)exp(tB) = (I+tA+t?4?/2+0(t?)) x (I tB - C B?/2 + O(€)) 
I +t(A+B)+t?(A?+ АВ + B?) + O(£). 


Collecting all the terms involving non-zero powers of t into a single term X we have 
exp(tA) exp(tB) = І + X. For t, and hence || X||, sufficiently small we can use the series 
expansion log I + X) = X — X?/2+ X3/3—--- to get 





C(t) = log(I +X) 


t(A + B) +t?(A? +2AB + B2)/2 — t?(A + В)?/2 + O(£). 





This expression clearly satisfies (L28) up to O(t?). The series for C(t) converges as it is 
obtained via elementary operations of convergent series. The goal of the Baker-Campbell- 
Hausdorff (BCH) formula is to find an explicit formula for the coefficients of the series for 
C(t), and to express the coefficients in terms of commutators acting on the matrices A and 


"More generally, one can define the exponential map in the same way when A is an element of a Lie 
algebra, for example, the linear operators over a vector space which need not be finite dimensional. 
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B. In order to do this we first give a few useful definitions and results. 


The Bernoulli numbers 


We will later find it necessary to make use of the Bernoulli numbers. We consider the 


function 
2 
НЕБЕ БЕЗ € C; 1.29 
exp(z) — 1 ? d 





the possible poles of which are z = +27in, n € №. Expanding (120) as a power series gives 


1 
oo 


2 “ү & 
== у= ——— = Юуг“ (1.30) 
Sea), а 
exp(z) – 1 = (7 +1)! -— 
where the Dj, are unique and of the form D; = B;,/k! with Bj € Q. The Вк are the Bernoulli 
numbers, the first few of which are given by Во = 1, В! = -i Ba = c B, = 3 Bg = 


b and Bj = 0 for k = 2n +1, n € N. Further values of Bẹ can be found using the recursion 


k—1 
—1 k+1 
B = — У By. 
, E» " 


The adjoint operator 


The Lie bracket [U,V] of two operators U and V is given by [U,V] = UV — vul If U and 
V are matrix operators then [U,V] is the matrix commutator and is zero if and only if U 
and V commute. If we take U to be fixed then for an arbitrary V of the same dimension 
the commutator [U,V] defines a linear operator V +> [U,V]. This is known as the adjoint 


operator and is denoted by ady(V) = [U,V]. We will use the notation adj, to denote 


*We will use the notation UV to mean the composition U o V in the case of U and V being infinite 
dimensional operators. When U and V are (square) finite dimensional matrix operators, (possibly resulting 
from the discretisation of an infinite dimensional operator), then UV simply denotes the matrix product. 
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repeated application of the adjoint operator, defined recursively via 





аа? (V) = V 
adi(V) = [U,V] 
аа? (У) = [U,[U,V]] 
ай (У) = [U,[U,...,[U,V]...]] 
i—times 
= [Шаа р 


The differential of the exponential map 


The differential of the exponential map is given by 


Sin : (1.31) 
ағ is 


We will show that this is equivalent to 


dexpy (V) exp(U) (1.32) 
where m 
1 
dexpy(V) = 2 трата (У): 


Expanding (L.31] in a series and then differentiating with respect to t and evaluating at 

















t = 0 gives 
— exp(U + tV) = —| 1+(U+tV)+—(U*+t(UV+VU) + O(t^)) 
dt = del; 2! 
1 
e (U? -t(U?V +UVU + VU?) + O(£)) +: 
1 
u (U" + t(U"7 V +U" °VU + + VU") + O(8?)) ++ 


1 1 
= Var (Ve VU eV ae OVE eV ee 


1 
des (U"-1v cU"? VU +... + VU") $e, (1.33) 
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Expanding (1.32) gives 


1 1 
dexpp(V)exp(U) = (1 + gradu + gradu Tee ) exp(U) 


V+ (UV - VU) + = (U[U, V] — [U, V]U) 4--- ) enn 


1 1 
= (v te UV = VU 5] (U?V -2UVU + VU?) +-- ) exp(U) 


(1.34) 


In order to simplify the task of matching the terms in (1.33) and (1.34) we use the notation 
introduced in [Tse, 2003]. Given that the matrices U and V need not commute and that we 
wish to compare terms involving the matrix V multiplied from the left and right by some 
power of the matrix U, we consider a term of the form U^—-!VU?. We let Z” denote a 
term with a total power of n for the matrices U and V. We let Y? signify that the same 
term is multiplied from the right by a j—th power of U. We also use the convention that 
Y? = І this gives, for example: 


y = zgy?-z 

UV = 22ү = 2? 
VU = ZY 
UVU = ZY 
Ue yi s Zen, 


It is easy to see that this notation gives a unique representation for each term in (1.33) and 
(32. Rewriting (1.33) in the new notation gives 


d 1 1 
a oP +10) = E IST) л с 
t=0 i 
1 n n nyzrn—1 
+ (2% + Z"Y 4 ----Z^y?" 71) 4... 
n! 
1 1 
= De Ur He PP e spen 
1 
a Нан 


In order to rewrite (L34) we note that the exp(U) term in the expression is multiplied to the 
right of the expansion for аехрг; (У) and hence we can rewrite exp(U) as exp(ZY). This 
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gives 
dexpy(V)exp(U) = (2+ 5 +-(2- 72ү) + Dd – 273ү + ZY? +... ) exp(ZY) 
= a" Z*(-Y)4 2 ZU - Ү)? + ees) 
= (1+zu-n Y)+ an yy +g- +) 
x(I- Y) Nd 


However, since the Taylor expansion of exp(ZY) contains only terms of the form (ZY)™ 


(i.e. only powers of U), it commutes with (7 — Y)~! and the last line above becomes 
(exp(Z) — exp(ZY)) ü - Y) 
We can then expand the exponentials in series to get 
(e? – е27) (1 Ү)-! = (1z tZ ) — (rave zvte )) (=): 


(21 = Yez OU = Y?)+5 ZU yet, ) (T-Y)! 


1 1 

Ze res 
1 

СРВ eee 


where the last equality follows from using a geometric series to write 


(Т — Ү")( 2» Yi. 
Hence, (1.33) and (34) are equal and 


= dexpy (V) exp(U). (1.35) 


d 
d exp(U + tV) = 
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The inverse of dexpy 


We also want to find an expression for dexp;;', the inverse of dexpy. To do so we observe 


that if A is an eigenvalue of the ady(V) then the eigenvalues of 


are given by 


Hence, for А Z +2rin, n € №, dexpy(V) is non-zero and has an inverse such that 





оо 


24 
dexpy (V) (> mn) S (1.36) 


k=0 


Comparing (1.36) with (1.30), from the definition of the Bernoulli numbers, it is clear that 


dexp; (V) = У —adı(V). (1.37) 


The BCH formula 


We now have the necessary tools available to give an explicit form for C(t) in (28), the 


BCH formula. We follow the proof given in |Hairer et al., 2002] and |Varadarajan, 1974| and 


show that if C(t) satisfies (L.28] it is also a solution of the initial value problem 


Zen =A+B+ А — B,C(t)] + > adi (А +В), C(0)=0 (1.38) 


We consider 
exp(sA) exp(tB) = exp(Z(s, t)) (1.39) 


with a smooth matrix function Z(s,t) and small values of s and t. Using (35) we take the 
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derivative of the left- and right-hand sides of (1.39) with respect to s. This gives 


Aexp(sA)exp(tB) — 2 exp(Z(s,t)) 
д 
= 4ехрды (3526.0) exotztos t) 
(inserting (L39) ) = dexpz(, (226) exp(sA) exp(tB) 


д 
> A = dexpz(s1) (26.0) 
we then use (L37) to find the inverse of dexpz(,, which gives 


д = c Br. 
752 (st) = dexp7/, »)(A) = Tr adze (4). (1.40) 
We now take the inverse of (L.39), exp(—tB) exp(—sA) = exp(-Z(s,t)), and proceed in a 


similar fashion, taking the derivative of both sides with respect to t 


Bexp(-tB)exp(-sA) = © exp(—Z(s,1) 
д 
= —dexp_7(s 1) 7260) exp(—tB) exp(—sA) 
д 
= В = dexp zu) (260) 


then applying the inverse of dexp_ хү, 





д Е By 
a (n = dexp (В) = 57 8d zo.) (B) 
=0 
1 X By 
= B — 8d z(s) (B) + У, a zen (B) 
k=2 
1 c Bi 
= В+ Zadzs(B) + qp adze) (B) (1.41) 
k=2 
Where the last equality follows from —3ad.z = —i(-ZB — B(-Z)) = $adz(B) and 
ad? z = (-1)*[Z,|...,[Z, B]...]] = (-1)*adE(B), however, В, = 0 fork =2n+1, neN 


hence Bad% (B) = ByadE(B) for k > 2. 


Comparing (£28) with (£39) we see that exp(C(t)) = exp(Z(s,t))|,_, and hence C(t) is 
the solution to ic) = £Z(s,t)|,_, + #Z(s,t)|,_,- It follows upon substitution of (L40) 
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and (LAT) that 
a = A—- = = a 
4:00 | 512, А] Ы a и) Ы |» + 2[2,В] > Ar sin) 
s= = t 








= k 
= АВЧ TA- B,C(t 1+ у; Жай (o (A +B) 
k=2 
hence C(t) is the solution of the initial value problem given by (38). 
We now express C(t) as a Taylor series; C(t) = XX] Cpt”. Differentiating with respect 


to t gives 


d оо 
—С(%) = C + 2С5® + 303 +- = kc, t*-}. 1.42 
ay © Mt) = Ci + 2Cat + 363 + Ук, (1.42) 


Substituting C(t) = 0%, Cpt" into the initial value problem (38) and comparing the 
coefficients with those of (1.42) gives 


Cı = A+B 
C; = 4[A-B,A+B] = 1А, В] 
Сз = $[A- B,SIA 8] = 34,4, B] + 75[B,B, 4] 





Cy = А,В, В, A] 
1.43 
Су = [A. A. A, А, B] – =;[В, В, В, B, А] ve) 
+ЫЛА,В,В,В, А] + 34518, А, A, A, B] 
+т[А, A, B, B, A] + Zu [B, B, A, A, B] 
where we have used the notation [A,..., B,C], to represent [A,...,[B,C]] and where the 


expressions for C4 and Cs are simplified using the Jacobi identity [A, |B, C]] + [B, [C, A]] + 
ІС, [A, B]] = 0. The matrix C(t) in the BCH formula (1.28) is therefore given by 


C(t) = (A+ B) SIA, BP + (А, A, B] + [B, B, A) O +14, B, B, Alt +O) (1.44) 


2: Splitting Methods for Nonlinear 
PDEs 


In this chapter we use the nonlinear Schrödinger equation (NLSE) to show how splitting 
methods can be applied to nonlinear oscillatory PDEs. A split-step method for the NLSE is 
developed with time propagation accomplished by first order Lie splitting while for the spatial 
discretisation we consider both a finite difference and a pseudo-spectral Fourier method. 
The Fourier method requires that periodic boundary conditions are used for the spatial 
domain and hence we investigate the properties of both versions of the splitting method 
applied to problems with a periodic solution. The methods are studied with respect to their 
treatment of conserved quantities, dispersion and instability. Step size restrictions are given 
which prevent numerical instability for the periodic solutions. In order to be able to better 
characterise the behaviour of the methods, the analytic properties of the NLSE are also 
studied. 


2.1 Splitting of nonlinear PDEs: Time discretisation 


We consider a nonlinear PDE of the form 
ш = би + N (u)u, (2.1) 


where £ is a linear operator, such as the Laplacian, and N (u) is a nonlinear term. Given 
a solution u(z,t4) to ŒI) at time tn we wish to find u(x,t, + т) the solution at the next 


time-level, where 7 is the (positive) temporal step size. To this end we use the following 
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formal solution to { Т) 


ЕТЕ "ET ( | Е NN as) ale) 


tn+T 
= exp (re + (u(x, sas) ubt tn). (2.2) 
tn 
Approximating the integral in 2.2] we get 
u(x,t, +7) = exp (r£ + TN (u(z,t,,))) ulz, tn) + О(т?). (2.3) 


In general (2.3) is first order accurate, however, under special circumstances it is exact. In 
particular, for the case of the NLSE where N(u(z,t)) = iq|u(z,t)|? we can evaluate the 


nonlinear term exactly. To see this we consider the nonlinear differential equation 

Ou | 2 

S —(a-ib)ulu, a,bER,  u-u(zt). (2.4) 
We then define y = |u|? = uu* which gives 


ду ди, Qu* —— ; 4 . 4 
Z = Sut tut = -(atib)lult- (a ib)lul 
= — 2ау°, 





so that &l = 2a. Since a is a constant we integrate with respect to t and get 


ol 1 
—-dt=— = [mats +c 
Oty y 
and thus 
1 1 
= — — = 2a(t — to), 
y Yo 





u(t = to)|?, hence y = yo/(2ayo(t — to) + 1) and 


where yo = [uo]? 





|u|? = [шо]? 
(2aluo|?(¢ — to) +1) 
This gives 
Ou | —(а + ib) |ug|? 
= = = b\jul? — ол. -., 
дї re 2а[ио|2(® — to) +1 


—0 
= -wrt _ 
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where 0 = (1 + ib/a) and y = 1/(2aluo|?). One can then verify that the solution to (25) 


and hence to (2.4) is given by 


1 1 
u(t) = uo exp (ошо) ехр (зове — to + ») А 
and, therefore, 


u(to + T) = uo exp (-5 (1 + iz) In (2a|uo|?r + D) . 


We then expand the In (2a|uo|?r + 1) term as a Taylor series to get 
и( +T) = Wexp (-3 (: + i) (satur — 5 alus c? Te )) 
= ugexp (—ib|uo|?r + O(a)). 
Taking the limit as a tends to zero we find 
lim u(to + т) = ug exp(—ib|uo|?7). 
From this we see that the solution to 
pu = iq|u(z, t) |?u(z, t) 


ot 


is given by 
u(z,ts +T) = ue, tn) exp(iq|u(z, tn)|?7). 


Since (2.6) also has the formal solution 
tn +T 
u(x, tn +7) = u(z, t4) exp (af [и(2, as) А 
tn 
we can compare (2.7) and (2.8) to see that 


in+T 
f ы р =н ДЫЙ» 
tn 


which in the case of the NLSE gives an exact expression for the integral in (2.2). 


(2.9) 


In the case of different nonlinear terms other approaches are available. Many of these 


involve approximating the integral by using a predictor-corrector type approach to first 


calculate an approximate solution at some intermediate point(s) and then using a quadrature 


rule to replace the integral. 
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To obtain a first order splitting method we need to replace the right-hand side of (2:3) 
with an expression which we are able to calculate explicitly. We do this by applying the 
BCH formula as given by equations (1.28) апа (44). From this we see that a first order 
approximation to the right-hand side of (2.3) is given by 


exp(r(£ +N (u))) = exp (T£) exp (TN (u)) + O(r?) (2.10) 


which is known as Lie splitting] 


In order to implement the time discretisation we define U™ (x) to be the discrete approx- 


imation to the solution u(x,t) at t= mr, m € №. We also define the intermediate solution 


И (а) = exp (TN (U™ (xD) U” (x). (2.11) 
Equations (2.3) and(2.10) then give 
U™H (x) = exp (r£) W” (£). (2.12) 


Although we refer to W(x) as the intermediate solution it is important to note that it is 
not an approximation to the theoretical solution but rather is the solution after solving for 
the effect of only the nonlinear operator N. We also note that if U(x) is L-periodic then so 
too is W™ (x). This is important if we wish to use periodic boundary conditions as it means 
the boundary conditions are not violated by the introduction of an intermediate step. 

It is easy to see that if the time step is imaginary then the scheme given by (BIT) and 
(2.12) is unitary. That is, if we define y(r) : C^ — С, (т) : U"(z) 5 U™ 1 (2) then 
y*(r)y(r) = I. This has consequences for the conservation properties of the numerical 
scheme as we will see later. 

The NLSE is an example of a PDE of the form { Т), with £ = iV? and N (u) = iq|up?, 
to which we can apply the splitting method above. First, however we look at some of the 


analytic properties of the equation. 


2.2 Analytic properties of the NLSE 
The one-dimensional cubic nonlinear Schrédinger equation (NLSE) is given by 
iu; + Ure + ии = 0 (2.13) 


‘Or Trotter splitting after [Trotter, 1959], or Kato splitting, or any combination of the three. 


2.2. Analytic properties of the NLSE 25 


where u = u(x,t) is a complex valued function of time t and space x defined for t > 0 and 

x € R and q is a real-valued coefficient. For convenience we rewrite equation (2.13) in the 
А 9? . . 

form of (ZI) by setting С = Zi and N (u) = iqlul?. 


2.2.1 Conserved quantities 


The NLSE satisfies an infinite number of conservation laws. The “mas of the solution is 


given by 


M(t) = ГА Mem ces (2.14) 


—oo 
If we consider u(x,t), a solution of the NLSE, with support on the interval (a,b) then 
calculating the derivative of (2.14) we find 


b 
T | (иш + ши*)ах 
b 
= / (u(—iux,, — iq|u[?u*) + (ius; — iq|u|^u)u*)dz 
a 
b 
= if (изаи — Ux, )аз 
a 


b b 
| b EN De 
= dugu'|, — | Uuzde — iiuzul, + if Ui 
a a 
= 0. 


We therefore write the mass conservation law as 


а T 2 

T ju(x,t)| dz = 0. (2.15) 
—оо 

In the case that the solution is an L-periodic function in space such that u(x-- L, t) = u(x,t), 

then the integral in equation (2.15) can be restricted to a single period of length L rather 


than the bi-infinite domain. 
Other physically important conserved quantities of the NLSE are the “energy” or Hamil- 
tonian given by 


H(t) = ГА (Ius (x, t)? _ |u(a, t)|*) dz, (2.16) 


—oo 


?'The same quantity is also often referred to as the energy of the solution. The term “particle number" is 
also occasionally used. 
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and the “momentum” or “current” given by 


Te / "utet te а Pete а, (2.17) 


00 


Here we will restrict ourselves to dealing with the conservation of mass given by equation 


(2.15). 


2.2.2 Dispersion 


The equation 
u(x,t) = Aexp(i(ka — wt)). (2.18) 


describes a periodic wavetrain with constant amplitude A, frequency w and wave number k, 
(k =1/A, А = wavelength). By substituting equation (2.18) into (2.13) it is easy to see that 
(2.18) can be a solution of (2.13) only when 


ш = k? — ql A}?. (2.19) 


For positive values of q this dispersion relation implies that at a fixed frequency the wave- 
length of the solution becomes shorter as the amplitude of the solution increases hence when 
q is positive N (u) is known as an attractive nonlinearity. The case for which q is negative 
and where the wavelength of the solution increases with increasing amplitude is known as 
a repulsive nonlinearity. The presence of the amplitude in the dispersion relation (2.19) 
is typical of nonlinear waves and can be viewed, for example, as modelling the nonlinear 
amplitude response of the refractive index of an optical fibre. As we will see in the next 
section, when q > 0 the dispersion relation can lead to analytic instability of solutions of the 
NLSE. 


2.2.3 Instability 


The NLSE can give rise to analytically unstable solutions. The interaction between this 
instability and the conservation laws, particularly the mass conservation leads to interesting 
behaviour such as “recurrence”. The instability of periodic solutions of the NLSE can be 


studied by considering a perturbed solution of the form 
u(z,t) = (1+ (2, t)) ü(x, t) (2.20) 


where je]? < 1 and 
ü(r,t) = Aexp(i(kx — wt)) (2.21) 
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with A,w and k chosen such that the dispersion relation (2.19) is satisfied. From (2.20) we 


calculate the following terms up to first order in e: 
juu = (1te+e*)|a|?u, |2 = |A]? 
= |u Paa +e) + |Al?a(e + e*) 
up = Ut + Eù + Ete 
Ure = Ung + 26505 + Exel + Ейр. 
Substituting these into equation (2.13) and simplifying using à, = ikü gives 


(1+6) + &ü = (їй + iqt t) (1+) + (—2Кє„ + ig| Al (e + €*)) à. (2.22) 


Since the coefficients A, k and w in à were chosen so as to satisfy equation (2.19), à is a 
solution of (2.13) and equation (2.22] reduces to 


Et = leg, — 2ke, + ig| Al? (e + e*). (2.23) 


If we assume that =(2, #) is periodic in space on the interval [- 52, 2L] then we can write = 


as a Fourier series, 


оо 
2тт› 


elx,t) = > En(t)exp(itnt); Um = En (2.24) 


n=— со 


It is then possible to use equations (2.23) and (2.24) to show that 











d 2 ê 
dr Kd = Gan NA , N= —OO,...,00; п #0 
E n = Tt. 
where 
g, —{| ФАР = 2kun = An, gl Al? 
= ; 
-q|AP, gl Al? — 2kun + u2 





The eigenvalues of Gn are given by 


An = —2їЁиһ + л 29 | Al? — u2 


and therefore Gn has an eigenvalue with a positive real part when 





0 < u2 < 2g| Al. (2.25) 
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Clearly, equation (2.25) is satisfied only when q > 0. In this case the vector [és el grows 
exponentially with time. The growth in the Fourier coefficients for which equation (2.25) 
is satisfied results in the growth of the corresponding frequencies in the perturbation and 
hence, leads to instability of the analytic solution. 

Although the analytic instability leads to growth of arbitrarily small perturbations when 
q > 0 and (2.25) is satisfied, the mass conservation of solutions of the NLSE means that this 
growth is held in check. This can lead to a phenomenom known as recurrence where unstable 
solutions first grow but then return to a state close to the initial conditions of the system as we 
shall see in section [2.5] This phenomenon was first reported in [Fermi et al., 1974]. An early 
example of a numerical study of recurrence in wave equations is given in |Zabusky and Kruskal, 1965 

For the case considered here where the potential term of the NLSE is given by A/(U"" (x)) 
= iq|U"(x)|? we can apply the previously mentioned splitting method and calculate the 
intermediate solution W™ (x) without any further approximations. However, the calculation 
for the final solution U™*+(a) involves the Laplacian operator С which must be discretised. 
One possibility for doing so involves using second order central differencing in space to obtain 


a finite difference approximation, we study this next. 


2.3 Spatial discretisation — Finite differencing 


Although the analytic solution of the NLSE is defined on the real line, any numerical treat- 
ment requires that we impose boundary conditions at the endpoints of some finite interval 
in space. In order to be able to compare the results of the split-step finite difference method 
with those of the split-step Fourier method, and because we initially wish to derive results 
for periodic solutions, we use periodic boundary conditions on the interval [— iL, iL]. It 
is, however, easy to apply the finite difference discretisation to problems with, for example, 
Dirichlet or von Neumann boundary conditions. 

We divide the interval [-iL, ŽL] into N equal subintervals with the grid spacing given 
by h = L/N and grid points z; = jh, j = —N/2,..., N/2. The approximation to the 
solution u(xj, тт) is denoted by Ui". We then use the (1,1) Padé approximation to the 


exponential function to get 
| in) L 
explirL) = |1— jc I+ git (2.26) 
where J denotes the identity орао By comparing the Padé approximation with the 


3 I will also be used to denote the identity matrix, the intended meaning should be obvious from the 
context of usage. 
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Taylor expansion of the exponential function it is possible to show that the leading error 


term in (2.26) is i7? £?/12. That is, the approximation is second order accurate. Setting 


(2.26) into (2.12) we get 
(1 Е size) pneu = (1 + ge) W(x). (2.27) 


We now discretise the Laplacian operator and the space variable by replacing £ with the 


second order central difference operator £; defined by 


1 2 —N N 
We also set 
UT. ,—7Ux , and Un —UX 
2 2 2 
to ensure that the solution is periodic. 
We can now write (2.27) in matrix notation; 
1. m+1 1. m 
I - girs U™™ a) =( i+ 575 W™ (x) (2.29) 
where 
wW = exp(irq|U7 |)U7" (2.30) 


with И" = D. Lu QU AI т = 7/h? and 
2 2 


The temporal discretisation given by (2.11) and (2.12) is of order O(7) while the spatial 
discretisation given by (2.28) is of order O(^?). The Padé approximation in (2.26) is accurate 


to second order in тС and hence has accuracy of О (°) The overall accuracy for the 


split-step finite difference scheme is therefore O (7 +h? + ()’). The r/h? term means 
that not just 7 and h but also the ratio between them is important in determining the 
overall error of the scheme. We will treat this in more detail in section [2.3.3] when we 


consider the instability of the split-step finite difference scheme. 
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Although the split-step finite difference scheme given by equation (2.29) is implicit, it 
is possible to make use of the quasi-tridiagonal structure of the finite difference matrix S 
and employ an optimised version of Gaussian elimination to solve for U™*+. In each step we 


have a quasi-tridiagonal system of the form 


di uy q gi b 
lj. da up T2 bə 
з. : = 
UN-1 ZN-1 bn-ı 
p In-ı dn XN bn 


Treating each element individually to avoid multiplying by zero we proceed by eliminating 
first the lower diagonal elements and then the p element in the lower left corner. This gives 


a system of the form 


|^ ui q | | X1 | by | 
| d? ul p* 29 by 
ED a TN-1 m 
di; EN bn 


where the crosses indicate values which have been altered. The solution can then be obtained 


via back-substitution. 


2.3.1 Mass conservation 


In order to show that a discrete analogue of the total mass is conserved by the split-step 


finite difference method we rewrite equation (2.29) as 
Um. уут = S(su™ + SW”). (2.31) 


We now pre-multiply this expression by (U™+! + W™)Ħ# and take the real part of both sides 
to get 


Re [Zom j wn) sie ж и") Ве (те d wt = w™)} 


(UHHH ymt = (Ш"УН г" (2.32) 
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where the second equality follows from using equation (2.30) which gives 
ww" = (exp(irglU™|?)U™) exp(irg|U" PU" = (07) 0". 


Since the finite difference matrix S is real and symmetric it has only real eigenvalues hence 
for any complex vector y we have Im{y Sy} = 0. It follows that the left-hand side of (2.32) 
vanishes and that (U™+!)#U™+! = (U™)#U™. The split-step finite difference scheme 


therefore satisfies a discrete version of the mass conservation law (2.15) given by 


N N 
Nj Nai 


h M qUPTP-h MS UP, тел. (2.33) 
ј= =“ j==X 


2.3.2 Dispersion 


In this section we will see that using the split-step finite difference scheme it is possible 
to obtain a discrete analogue of the nonlinear dispersion relation given by equation (2.19). 
As in section [2.2.2 we assume a periodic solution of the form Ui" = Aexp(i(kx; — wmr)). 
Substituting this into the right-hand side of equation (2.29) and using (2:30) gives 


(1 + Zs) И" = (1 + Zs) exp(ikz) A exp(irq| A|?) exp(—iwmr) 
=R 
N_1 
where 2 = le; у: Since x; = [h we have 
25 5 
Г = exp(iklh) + 5 (exp(ik(! — 1)h) — 2exp(iklh) + exp(ik(l + 1)h)) 
= exp(iklh)(1 + ir(cos(kh) — 1)). 


Setting s = cos(kh) — 1 we get 


(r^ 25) we] = RjAexp(irq|A|?) exp(—iwmr) 
l 


exp(iklh)(1 + irs) Aexp(irq|A|?) exp(—iwmr). 
Similarly, [(7 — #5) Un] = exp(iklh)(1 — irs) Aexp(-iw(m + 1)7). Equation (229) now 
gives (1 — irs) exp(-iwr) = (1 + irs) exp(irq| A|?) and hence 


_ 1—78 
С L+irs’ 





exp(ir(w + а|А|?)) (2.34) 
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Setting T(w + q|A|?) = —20 we rewrite equation (2.34) as 


cos(—20) + isin(-20) = ————5 — 


Taking the real part then gives 


1— rs? 


= - 2 — 2 
cos(20) = 1 — 2sin^(0) = 2 cos“ (0) – 1 = Pape 


from which it follows that sin(0) = rs/v 1 + r?s?, cos(0) = 1/V1+r?s? and 


ш = = — 4А. (2.35) 
In order to compare the dispersion relation given by (2.35) with that given by (219) 
we take the limits h > 0, т — 0 such that r = r/h? remains constant. This gives 
s = cos(kh) — 1 = —k?h? /2 + k*h*/24 — --- and since s? — 0 as h — 0 we have авт — 0, 
0 ~ sin(0) = rs/V1+r?s? c rs and therefore 


1 
0 = = + O(k*n*) 
where the constant implied by the O(k*h*) term depends on r but not on h or т. Hence, 


ш = zer + O(k*h5)/7 — А> = k? — q| Al? + O(k4h?) (2.36) 
and the dispersion relation given by equation (2.35) tends to the analytic dispersion relation 
(2.19) as h and 7 tend to zero. From this we see that the accuracy of the discrete dispersion 
relation depends on whether kth? is sufficiently small. Recalling that k is the wavenumber of 
the solution it is clear that this requirement is none other than the intuitive notion that highly 
oscillatory problems should be solved on a fine spatial mesh. The approximation is exact 
when k — 0 and hence in the absence of rounding errors the finite difference discretisation 


is exact for solutions independent of =. 


2.3.3 Instability 


Here we are interested in the ability of the split-step finite difference scheme to model the 
analytic instability of solutions as discussed in section [2.2.3] and whether new instabilities 


without analytic counterparts are introduced. We consider a discrete solution U” which is 
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independent of space and given by 
U” = aexpl(ig| Amr). (2.37) 


This corresponds to the solution obtained from the temporal discretisation of (2.21) with 
k = 0 and w chosen so that the dispersion relation (2.36) is satisfied. We now perturb (2.37) 
to 

UP =U™(1 + ef") (2.38) 


where lew" |? < 1. Setting (2.38) into (2.30) and retaining only terms up to first order in &7' 


gives 


WP = explirgÜP(L+EP))UR(1+ ER) 
= exp(irg|Al?(1 + Er + €?*))UT (1 + e?) 
= exp(irq|A|) exp(irq(e5. + ej" ))A exp(itq|A|?m)(1 + 7) 
= UMM +irg(e” + e7*) + O((7)?))Q + =") 
~ О e + irge? + eT"). (2.39) 


Substituting (2.39) and (2.38) at t = (m + 1)т into (2.29) then gives 
(1 = Sch) gH = (1 ү 244) ((1 + iq] Ar)? + igl A]? rer) (2.40) 


where Lp is the second order central difference operator defined by (2.28). Assuming that 
€" is band-limited and periodic on the interval [-L/2, L/2] allows us to write it as a discrete 


Fourier series; 
X- 
2 
2тт 


£j = > En exp(iunz;), Hn = T (2.41) 
N 


Nn=- 





Setting the Fourier series (2.11) into (2-40) and collecting terms in exp(+pnx;) gives 





А m+1 " m 
gt. Ё n 
where 
ca d, (1 + ig|A|?7) d,iq|A[^r f= 1+ rsy 
i -diigAPr  d;(l-igAPr) |' ^" 1—з„' 








Sn = cos(us h) — 1 and r = т/ћ?. 


34 Chapter 2. Splitting Methods for Nonlinear PDEs 


The eigenvalues of Gn are given by 





An = bn + \/82 —1 (2.42) 


where 


_ .1- r?s? – 2q| Ars, 
Е 14+ 7282 | 


Bn 


ak 


We see from this that the Fourier modes [2,,,2* „| will grow exponentially if one of the 


eigenvalues given by (2.42) has |А„| > 1, ie. if |8,| > 1. Since s, < 1, this occurs when 
—rs—2q|A|?r < rs < 0 which, using the identity cos(z) — 1 = —2sin?(x/2), can be rewritten 
as 5 А 

ql Al? > ie) >0 (2.43) 


and similarly, when —2q|A Ав. < —2 which we can rewrite as 


2g APT, sin? (an) <-1. (2.44) 

Condition (2.13) can only be fulfilled when q > 0 and therefore corresponds to the 
analytic instability of section [2.2.3] For small Р, condition (2.43) gives 0 < и? < 2q|A|? + 
O(h?) which matches the condition for instability of the analytic solution, (2.25). The 
number of unstable modes in the Fourier series for the perturbation = is the largest whole 
number n such that condition (2.43) is fulfilled and, hence, all modes are unstable if n = N/2. 


'This situation is characterised by 


dle 


h> (zi) . (2.45) 


If we consider a small random perturbation, such as a rounding error, then for А large іп 
the sense that condition (2.45) is satisfied all Fourier modes of the perturbation may grow 


exponentially and the numerical solution will display severe instability. 


The second condition for instability given by (2.44) is only possible for q < 0 and hence 
is a purely numerical phenomenon. Such instability can be prevented by always choosing т 


such that 
h 


< —- 
v2lallAl? 


which ensures that none of the Fourier modes are unstable. 


(2.46) 


While we have seen that |3,,| > 1 leads to instabilities, for |G,| < 1 we have Bn —1 <0 
and equation (2.42) gives |A|? = 82 — (82 — 1) = 1. That is, the method is neutrally stable 


2.4. Spatial discretisation — A Fourier method 35 


when neither of the instability conditions are fulfilled. 


2.4 Spatial discretisation — A Fourier method 


If we have periodic boundary conditions it is possible to replace the finite difference discreti- 
sation of section [2.3] with a Fourier method. This method has the advantage that calculating 
the Laplacian in Fourier space introduces no additional error to the spatial discretisation and 
also that the necessary DFT/IDFT operations can be computed in only O(N log5(N)) oper- 
ations using a fast Fourier transform (FFT). Further computational gains can be made by us- 
ing a parallel implementation of the FFT algorithm (Briggs and Henson, 1995}/Van Loan, 1992]. 

To implement the spatial discretisation of the Laplacian using Fourier transforms we first 
replace U"**1(z) and W""(z) in equation (2.12) with their Fourier series. This gives 


amt = ехр(-йи2т) б", m= -00,...,00, aL 


where ü+! and i)?" are the Fourier coefficients of the spatially continuous functions U”+!(x) 


and W™(x) as given by (LIJ). We then replace 10" with its discrete equivalent W/" by 
calculating the DFT 


~ h : Е 
з= А E opima WP, n= By BA (2.48) 


where W? = exp(irq|U7" |)U7. We now advance the solution one time-step ahead in 
the frequency domain by calculating the discrete quantity UM*+! such that equation (2.47]) 
becomes 

Omt - exp(—ip27)W7", n= M. —1. (2.49) 


Finally, we return the solution to the spatial domain by applying the IDFT to (2.49) 


Nj 
2 
pre У ехр(іньх;) "і, j=... E - 1. (2.50) 
aN 


п= —5 


From equation (1.25) we see that if the function U™ (x) is band-limited with à7" = 0 for 
|n| > m,m € N then the Fourier discretisation of the Laplacian is exact at the points x; 


whenever N > 2m. That is, whenever the number of points in the Fourier discretisation is 


4 That is, no error other than that due to truncating the Fourier series by using a sum over a finite number 
of terms when calculating the DFT. 
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greater than twice the bandwidth of the solution. We will later see that in practice it is also 
possible to obtain good results using the Fourier discretisation method with non-periodic 
solutions. For example, in the case of a solution given by a single soliton it is possible to 
place the boundaries of the computational domain sufficiently far away from the peak of the 
soliton so that the solution near the edges of the domain is close to zero and hence periodic 


boundary conditions do not adversely affect the solution. 


2.4.1 Mass conservation 


We have already seen in section P.I] that the time discretisation given by (2.7 and (2.12) 
conserves the mass of the solution exactly. It therefore remains only to show the mass 
conservation properties of the Fourier discretisation of the Laplacian. This is simply achieved 


with the use of the discrete form of Parseval’s equality as given by (227. 
Using (2.48) and (2.30) with Parseval’s equality gives 


h = 
Уш? = унт? = Ay eee 
j j j 
We then use equation (2.49) and a second application of Parseval’s equality to show that 
h = " 
ТУЙ? = уйур = ушур 
J 1 J 


and hence, 


h^» qu'en pee (2.51) 
J J 


where j = —N/2,...,N/2—1 in all cases. Equation (2.51) is a discrete analogue of equation 
(2.15) and therefore, the split-step Fourier method exactly conserves a discrete analogue of 


the total mass. 


2.4.2 Dispersion 


Once again we assume a solution of the form 


U; = Aexp(i(kv; — wmr)) (2.52) 
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which implies 27/k-periodicity of the solution U7". Substituting this into equations (2.48) 


and (2.49) gives 


5—1 
wm = ; >, exp(-inaz;) exp(ir|UP PUP 

dec 
^2 
&—1 

= x У) ew — nz; + kag — wm) 
$n 
zer 

= (Ж (G(g|A?T — wmr)) n=1 

and 
Ümn- | 0 i 7 : 

Aexp(ir(g|AP? — om —k?)) n=1. 


Taking the inverse DFT of (2.54) then gives 


ao 


п Y Ur! = exp(i i(kz; + g|A|?r — отт — k?r)). 


-N 


coi 


The solution given by (2.52) evaluated at t = (m + 1)r is given by 
Ur = Aexp(i(kz; — w(m + 1)т)). 
This agrees with (2.55) in the case that 


w = k? — q] A}?. 


(2.53) 


(2.54) 


(2.55) 


(2.56) 


The dispersion relation (2.56) for the split-step Fourier method is therefore identical to the 


analytic dispersion relation (2.19). This is a direct consequence of the fact that the spatial 


discretisation is exact for N > 2m. Since the solution given by (2.52) has only a single 


frequency we have m = 1 and hence the dispersion relation is exact at the Fourier nodes for 
N > 2. Furthermore, since the solution (2.52) has |07? = |A|? the potential term N (U) is 


independent of of both x and t and the operators N and £ commute. The time discretisation 


given by (2.12) is also exact and hence the split-step Fourier method is exact for solutions 


of the form (2.52) in the absence of any computational errors. However, as such errors are 


inevitable we now study the stability of the split-step Fourier scheme so as to determine 
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under which conditions small perturbations can be magnified and lead to instability. 


2.4.3 Instability 


We follow the procedure of section [2.3.3] and consider a solution at time t = mr and 
independent of x given by U™ = aexp(iq|A|?mr). This solution is then perturbed to 
Up = um (1 + a) lew |? < 1. Again, we assume that the perturbation is L-periodic 
in space and express it as a Fourier series which we substitute into the expression for Wj” 
given by equation (2.30) and expand, retaining only terms up to first order in £7 to obtain 
the linearised intermediate solution Wj" = п +EP -Е?тд(є7' +e%*)). Substituting this 


into equation (2.48) then gives 








N-i 
> hix 2 . " ‚2тт 
И" = L X pet a Єў + ig| A|? T (e7 ae exp(-i725) 
i= 
_ f Or + (1+ ig AP rep + iglAPrer*)), n=0, (2.57) 
Un ( + ig APr)eT + ig|Al?re™)), n#0. | 
Calculating U7"*! from (2.57) via (2.49) gives 
рты Ј Un + (1 + tal A PTE + ig AP rep), n=0, 
" Ur (o ig| A[?7) exp(-innr)en + ig|Al?r exp(—iuzT)é7;)), n #0. 
(2.58) 
From equations (2.40]) and (2.41) we have 
[mal am+1 = 
бт UMM +E), n=0, (2.59) 
Отт, n 5 0. 


Comparing (2.58) with (259) we see that 2! = (1 + ig|Al?r)2” + 14| Al?rER”* and 


m 








where 


с | (Ф+їФАРт)ехр(—шут) iglAl’rexp(-innr) | 


—їч|А|техрбиҗт)  (1-— Ат) exp(ip7) 
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The eigenvalues of G, are given by An = Bn + 4/82 — 1 where 





Ba = (1+ %g|Al?7) exp(—ipnr) + (1 — ig|Al?7) exp(iunT) 


соз(илт) + g|Al?7 sin(unr) 
hence, instability occurs when 

| cos(u27) + q|A|?7 sin(u27)| > 1 (2.60) 
Multiplying equation (2.60) by 1/4/1 + q?|A|47? and setting 


1- Ф|Ат? 
20) = — 
сов(20) + par? 
(then using the identity cos(2x) = 1 — 2sin?(z) = 2cos?(x) — 1) gives 


1 q| А?т 


cos(0) = -h ae and sin(@) = -R A 
1+ q?|Al4r V1 + @ Alter? 


and hence, | cos(u27) cos(0) + sin(u27) sin(0)| > cos(0). Applying a second identity, cos(x — 


(2.61) 


у) = cos(x) cos(y) + sin(x) sin(y), finally gives 
| cos(u27 — 0)| > cos(0) (2.62) 
as the condition for instability of the split-step Fourier method. 
For the case of 3, > 1, the Taylor expansion of 3, as a function of 7 gives 
0 < u? < 24|А|? + Or?) (2.63) 


which, for small 7, is consistent with equation (2.25), the condition for analytic instability 
and, which holds only for g > 0. It can also be seen that this instability will affect the 
lower frequencies most strongly when it is present. The fact that (2.63) is independent of 
h, the step size of the spatial discretisation, is a consequence of the assumption that the 


perturbation &"* is L-periodic and therefore band-limited with m = N/2. 


In addition to the low-frequency instability given by (2.63), equation (2.62) allows for 
instability which occurs predominantly in the high-frequency modes. To see this we fix 
q|a?|, т and hence 0. We then define g(u?) = cos(u?r — 0) where u is treated as a continuous 


variable although in practise it takes only the discrete values given by un = 27n/L. Equation 
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(2:61) shows that 


de (0, 2), а> 0, (2.64) 
(350). q<0. 


The stationary points of g(u?) are given by u? = и» = (0--px)/r, p € No with the condition 
(2.64) requiring that р = 0 be excluded for the case of q < 0. It follows that |g(u2)| > cos(0) 
is equivalent to |u27 — 0 + рт| < |0| and therefore 


|0 


|н»— ир] << (2.65) 


the result being that whenever un satisfies (2.65), for some p, the n-th Fourier mode becomes 
unstable. For q > 0 and р = 0 condition (2.65) reads ти2 < |0| and we have instability in the 
lower frequencies corresponding to the analytic instability. The high-frequency instability 
resulting from p > 0, however, has no analytic counterpart and we therefore wish to prevent 


it. Since the largest value of u2 is given by 


we take р = 1 such that a sufficient condition for high-frequency stability for p > 1 is given 


by 
7? Ө+т 


h? T 


> IF (2.66) 


T 











When д > 0 we observe that for high-frequency instability 2 — ш? < 0 and so, making use of 
(2.64), condition (2.66) becomes rr?/h? — (0 + 7) < rn?/h?—- т < 0 and the high-frequency 
instability is prevented whenever 

h? 


m 2.67 
т< = ( ) 


Similarly, when q < 0 and —7/2 < 0 < 0, condition (2.66) implies that rn?/h?- (0+7) < 0, 


or rather, 





20m m oe 


ensures no high-frequency instability will occur. 


Although both (2.67) and (0.68) are sufficient conditions for avoiding numerical insta- 
bility, they are not necessary. For example, if |0|/7 is very small then it is possible that any 
u2 which satisfy condition (2.65) will fall outside the range of frequencies in the calculation, 


that is, |n| > N/2. In this case no numerical instability will be observed. 
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2.5 Numerical results 


In this section we show the validity of the stability conditions for the split-step finite dif- 
ference and Fourier methods which we presented in sections [2.3.3] and [2.4.3] for periodic 
solutions. 

For any method to provide qualitatively correct approximations to the solutions of the 
NLSE it is important that the method is able to reproduce the analytic instability of solutions 
and that numerical instabilities without an analytic counterpart are not introduced. In 
order to study these two properties with respect to periodic solutions for the split-step 
finite difference and Fourier methods we consider a solution of the same form as (2.18) with 
a = 1/2, k = 0 and w such that the dispersion relation (2.19) is satisfied. That is, a solution 


independent of space and given by 


TA Zep (=) . (2.69) 


With initial condition 1 
и(2,0) = 9 (2.70) 


Since the solution (2.69) is independent of space, both the split-step finite difference and 
Fourier methods are theoretically able to solve the NLSE exactly when the initial condition 
(2.70) is used. The theoretical bounds for the stability and instability of the solution given 
in sections 2.3.3] and [2.4.3] are studied by perturbing the initial condition (2.70) to u(z,0) = 
(1 + e(x,0)). We make use of the following three initial perturbations: 


1. 
E(x, 0) = 0. (2.71) 
In this case it is assumed that a perturbation will be introduced by rounding errors. 
2. 
e(z,0) = 0.1 cos = = <т< 5. (2.72) 


This perturbation introduces energy into only the fundamental frequency and will 
therefore be used to illustrate the low-frequency analytic instability predicted by equa- 
tion (2.25). The resulting initial condition looks like: 


0.6 T T T T 
0.4 1 1 А І 


-4 2 0 2 4 
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(2.73) 


d le), 0 


S gos 
0.1(1+ 22), #<r< 


This perturbation introduces energy into all the frequencies present іп the computation; 


the resulting initial condition looks like: 


0.6 ee o 
0.5 


0.4 - А 1 ‚ ‚ ‚ 4 














Because our interest here is in detecting when instability does and does not occur it does 
not make sense to compare the numerical solutions with an exact solution even though we 
have an analytic solution for the case when there are no perturbations. Instead we plot the 
modulus of the numerical solution and the modulus of the DFT of the solution. Since the 
DFT is symmetric, with |U^| =|U™,|, n = —N/2,..., N/2 — 1, we show the spectra for 


positive n only. 


2.5.1 Instability of the finite difference method 


The instability condition (2.43) corresponds to the case of the low-frequency analytic insta- 
bility. In order to verify this we use the perturbation (2.72) with the parameters q = 2, 
T = 0.005, L = 16 and N = 50. Condition (2.43) then becomes 


2x50? . of 2an 1 
0 < sin TARG 








162 2х 50 


which implies that the Fourier modes corresponding to |n| = 1,2 are unstable. 
Figure [2.T] shows that initially the |n| = 1 mode of the solution grows until t ~ 5 when 
energy starts to be transferred to the |n| = 2 mode. By t ~ 10 the solution has developed two 
distinct humps. The magnitude of these humps then decreases until at t > 18 the solution 
has returned to approximately the initial conditions. We see from this that the split-step 


finite difference method reproduces the analytic instability of solutions and gives rise to the 


well know recurrence behaviour documented in, for example, |Zabusky and Kruskal, 1965 
Fermi et al., 1974| [Yuen and Ferguson, 1978]. 


We also wish to investigate whether is is possible for similar instability to be introduced 
purely by rounding errors. To do so we take perturbation (2.71) and increase the temporal 
step size to т = 0.02. The resulting solution is shown in Figure [2.2] In this case no instability 


occurs until after t = 50. When the instability does occur it is the |n| = 2 mode which is 
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Figure 2.1: Split-step finite difference solution of the NLSE with perturbation (2.72) showing analytic 


low-frequency instability. A — 0.5, 


4= 2.0, І = 16, 


solution is shown above, the Fourier spectrum below. 


т = 0.005, N = 50. The modulus of the 
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100 





10 100 


50 


Figure 2.2: Split-step finite difference solution of the NLSE with perturbation { 717). А = 0.5, а= 
2.0, L=16, т = 0.02, N = 50. The instability is analytic and is due to rounding errors. The 
modulus of the solution is shown above, the Fourier spectrum below. 
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first affected before the solution returns to the approximate initial conditions after which 
the |n| = 1 mode become unstable. Because the instabilities in this case are introduced by 
rounding errors, the time at which they first occur and the order in which the modes are 
excited could depend on the method used for the solution of the linear system of equations 
given by 2.29). 

It remains to show that for the case of q < 0 condition (2.46) is sufficient to prevent the 
introduction of high-frequency instability with no analytic counterpart. To do so we take 
а= —8, L = 1, N = 20 such that condition (2.46) gives 

T< Su ER = 0.025. 
20,/2- 8/4 
We therefore consider two cases, one with т = 0.0245 for which we expect no instability, 
the other with r = 0.0255. For both cases we let rounding error be the only perturbation. 
For the case of т = 0.0245 we found that by t = 200 the change in the solution from the 
initial conditions was in the order of 107? and hence we conclude that any rounding errors 
introduced into the solution are not magnified. No plots are given for this case on account 
of their rather uninteresting appearance. 

The solution for the case of т = 0.0255 is given in Figure [2.3] where we see that high- 
frequency instability is introduced to the solution at t ~ 40. This is followed by a return 
to the initial conditions at t ~ 50 before the high-frequency instability reappears. The 
spectrum corresponding to the solution shows that this instability is nearly entirely limited 
to the modes with |n| = 9,10. This matches the prediction of equation (2:44) which for 
the parameters used here reads —4sin?(rn/20) < —1 and is satisfied for |n| = 9,10. The 
spectrum also shows that the instability remains limited to the high-frequency modes with 
no significant leakage into lower frequencies. This behaviour remained unchanged for longer 
integration periods with no notable increase in the instability of the modes |n| « 9 for 
integration periods in the order of 100,000 time steps, ie. t ~ 2550. We conclude from 
these last two experiments that the step size bound given by condition (2.46] is effective in 


preventing numerical instability for plane wave solutions with q < 0. 


2.5.2 Instability of the Fourier method 


We аге now interested in showing the accuracy of the bounds given by conditions (2.67) and 
(2.68). These conditions prevent the introduction of high-frequency instability for the cases of 
q > д апад < 0 respectively. To verify these conditions along with the condition for analytic 
instability given by (2.63) we use the perturbation given by (2.73) which introduces energy to 


all frequencies in the calculation and can therefore give rise to both high- and low-frequency 
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Figure 2.3: Split-step finite difference solution of the NLSE with perturbation (2.77 and A = 
0.5, q=-8.0, L=1, т = 0.0255, N = 20. The instability shown is numerical and is induced by 
rounding errors. The modulus of the solution is shown above, the Fourier spectrum below. 
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instability when conditions (2.63) and (2.67) are satisfied simultaneously. Choosing the 
parameters A = 4, q = 2, L = 16 and N = 32 condition (2.67) becomes т < 4 = 0.0796. 
If 7 is chosen near the limit of this inequality condition (2.63) predicts that low-frequency 








instability will occur in the modes corresponding to |n| = 1,2. We therefore consider two 
cases, the first with 7 = 0.0790 for which we expect only low-frequency instability and the 
second, т = 0.080, for which both low- and high-frequency instabilities can theoretically 
develop. 

This behaviour is confirmed in Figures 2.4] and 2.5] Figure [2.4] shows the modulus of 
the solution and its associated Fourier spectrum for the case of т = 0.0790. We see here 
similar behaviour to that shown in Figure [2.1] with instability occurring first in the |n| = 1 
mode followed by the |n| = 2 mode before the solution returns to approximately the initial 
condition at t ~ 26. This behaviour is then repeated with further recurrence of initial 
conditions at t ~ 46 and t ~ 66. It is clear, however, that in agreement with condition 
(2.67), no high-frequency instability is introduced. The fact that Figures [2.4] and [2.T] do not 
match exactly for the order in which modes become unstable and in the periods between 
recurrence of initial conditions is to be expected as in both cases the focus has been on the 
stability rather than the accuracy of the solutions. Choosing smaller values of 7 and h gives 


better agreement between the results of the two methods for longer integration periods. 


The solution for the case of т = 0.080 is shown in Figure [2.5] Until t > 40 the so- 
lution remains similar to that in Figure 2.4] with instability occurring in first the |n| = 1 
and then in the |n| = 2 mode. However, instead of displaying a second recurrence of the 
initial conditions the solution for 7 — 0.080 develops high-frequency instability in the mode 
[n| = 16. This is predicted by condition (2.62) which, for the parameters used here, reads 
| cos(7/8?n?0.080 — 0.0400)| > 1/4/1 + 0.0802/4 and is satisfied only for |n| = 16 when n 
is in the range of frequency modes present in the calculation. As in the previous cases, the 
growth of the unstable solutions, for both low- and high-frequency, remains bounded due to 
the conservation properties of the numerical scheme. 

We next consider the case of q < 0 for which the relevant stability condition is (2.68) 
and for which no analytic instability is expected. After setting q = —2, A = i, N = 32 
and т = 0.050, condition (2.68) predicts that for some L > Г, = 12.7848 the solution of 
the Fourier method will experience high-frequency instability while for L < L, the solution 
will remain stable. Further, condition (2.62) predicts that the instability will be limited 
to the mode n = 16. This predicted behaviour is verified in Figures [2.6] and [2.7] the first 
of which shows that for the above parameters with L = 12.6 and perturbation (2.73), the 


perturbation does not grow with time. 


Figure [2.7] shows the solution for L = 12.7. Until t ~ 24 the behaviour of the solution 
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Figure 2.4: Split-step Fourier solution of the NLSE with perturbation (2.73) showing analytic 
instability. A= 0.5, q=2.0, L=16, т = 0.0790, N = 32. The modulus of the solution is shown 
above, the Fourier spectrum below. 
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Figure 2.5: Split-step Fourier solution of the NLSE with perturbation (2.73) and A = 0.5, q = 
2.0, L=16, т = 0.0800, N = 32. The low-frequency instability is analytic, the high-frequency, 
numerical. The modulus of the solution is shown above, the Fourier spectrum below. 
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Figure 2.6: Split-step Fourier solution of the NLSE with perturbation(2.73) showing low-frequency 


numerical instability. А = 0.5, q = —2.0, L = 12.6, т = 0.050, N = 32. The modulus of the 
solution is shown above, the Fourier spectrum below. 


2.5. Numerical results 51 

















- 7, 


С 












AAS 
SN 


SS NN N 7 um IM 
= ROK a ~ | | 
m nn _ | \' I — ' 
IN c N cl hw 
SN AW WS 
SS NN | 


N 
SS 
Cos 





Figure 2.7: Split-step Fourier solution of the NLSE with perturbation (2.73) showing high-frequency 
numerical instability. А = 0.5, q = —2.0, L = 12.7, т = 0.050, N = 32. The modulus of the 
solution is shown above, the Fourier spectrum below. 
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is reminiscent of that for the stable solution in Figure 2.6] However for t > 24 the solution 
displays instability in the high-frequency modes. The Fourier spectrum confirms that this 
instability is indeed limited to the mode |n| = 16. 

These examples illustrate the accuracy of the conditions given in sections [2.2.3] [2.3.3] 
and [2.4.3] for determining when analytic or numerical instability will occur for plane wave 
solutions and show the effectiveness of the bounds determined in those sections for preventing 
the onset of numerical instability. They also show the effect of the conservation properties of 
the NLSE, with unstable solutions remaining bounded due to the exact mass conservation 


by both the finite difference and the Fourier method. 


3: Soliton Solutions of the NLSE 


Until this point we have only considered solutions of the NLSE which have been periodic in 
space. We now consider non-periodic solutions, in particular, soliton solutions. These are 


used to model a large number of physical phenomena from instability of deep-water waves 


to optical fibre commu 
nications and non-resonant plasmas 
[Maccari, 1999]. As well as being physically significant, solitons are important in their own 
right as a mathematical tool for the analysis of nonlinear systems as indicated by this state- 


ment by A. C. Newell [Newell and Moloney, 1992| who has published extensively on the role 


of solitons in mathematics, as well as their importance in nonlinear optics. 


“The soliton is to nonlinear science what the Fourier mode is to linear science, 


namely a fundamental ‘normal’ mode of propagation of a nonlinear system." 


Although the stability analysis of the previous sections focused on periodic solutions, we saw 
that for A — 0 there is neither analytic nor numerical instability, hence, for non-periodic 
solutions of small amplitude we can still expect stable solutions irregardless of the temporal 
step size. In this section we will restrict our experiments to using the split-step Fourier 
method. This may seem counterintuitive as the Fourier method imposes periodic boundary 
conditions on all solutions and can therefore affect the accuracy of non-periodic solutions. 
We will see, however, that by choosing a spatial domain wide enough so that the solution is 
zero in the vicinity of the boundaries we can still achieve accurate results. Additionally, the 
low computational cost of the Fourier method due to the highly optimised FFT algorithms 
available means that the method is significantly faster than the finite difference method even 
when compared with the specially optimised method for the system of linear equations which 
was outlined in section [2.3] For example, the calculation for Figure B.T]took roughly half the 
time to compute with the Fourier method compared to the finite difference method. Similar 


1 About 3N log,(N) complex operations for each step involving the linear term and for N = 2" spatial 
grid points. 
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results for the computational time required by the split-step Fourier method in contrast to 


a variety of other numerical schemes for the NLSE are given in [Таһа and Ablowitz, 1984|. 
An exact solution of the NLSE is given by 


u(x,t) = 2v exp(i(2¢a — 4(¢? — v?)t + (фо + /2)))sech(2vx — 8Gvt — xo) (3.1) 


where zo, v, ¢ and do are real constants relating to the initial position, amplitude, velocity 
and phase of the soliton respectively [Taha and Ablowitz, 1984]. We choose the constants to 
be то = 0, v = 0.5, С = 1 and фо = 0 so that the initial condition obtained from evaluating 
equation (3.1) at t = 0 is a soliton with amplitude one centred at x = 0. The result of using 
these initial conditions with the split-step Fourier method with N = 128, q = 2, т = 0.02 
and L = 50 for an integration period of T = 4 is shown in Figure [8.1] 

We see that there is good qualitative agreement between the numerical and the exact 
solution. The plot of the difference between the two solutions shows us that, as would be 
expected, errors are largest where the magnitude of the solution is greatest. This error is in 
the order of 107%, an observation which is confirmed by calculating the infinity norm of the 
solution over all иш which gives 9.87 x 107%. 


3.1 Moving computational domains 


In Figure [3.1] it is possible to see errors introduced by the periodic boundary conditions 
which take the form of a hump appearing near the left-hand boundary of the computational 
domain as the soliton approaches the right-hand boundary. This highlights the importance 
of ensuring that the domain is chosen so that the solution remains approximately zero in 
the vicinity of the boundaries. This can lead to the need for very large spatial domains and 
hence a large number of grid points when the integration period is long. 'To overcome this, 
it is possible to use a moving computational domain. 

'This involves shifting the spatial domain every few time steps so that the peak of the 
solution remains in the centre of the spatial domain [Knapp, 1995]. To do so, the values which 
correspond to points near the side of the domain furtherest from the peak of the solution 
(in the vector which describes the solution at time tn) are switched to the opposite end of 
the solution vector. Doing so ensures that the solution always has maximum amplitude in 
the middle of the solution vector. After moving the points in the solution vector it is also 
necessary to alter the values in the vector which corresponds to the grid points of the spatial 


domain on which the solution is defined. This is easily done by adding the number of points 


?That is, ||U — “||. = max; тахь |Uz" — u(hn, тт) 
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Figure 3.1: The magnitude of the numerical (left) and exact (right) solutions for the split-step 
Fourier method with initial conditions u(x,0) = exp(i(2x + m/2)))sech(z) and N = 128, q = 2, 
т = 0.02, L = 50, T = 4. The lower plot shows the magnitude of the difference between the two 
solutions over the same domain. 
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moved in the solution vector times the spatial step size to each term in the existing vector 
for the spatial domain. This ensures that the solution and the domain on which it is defined 


match at each time step. 


Given и, (solution vector at t = tn of length N), £n (space domain at tn), h (the grid 
spacing on the domain zn) and Isa (index such that |un(Iora)| = ||Un|loo) a sketch of a 
possible algorithm would be: Integrate for one time step to find ил+1 and find Inew such 


that |Un41Unew)| = ||Un+illoo) then, 


aI=I_new-I_old; 
if dI <= -1 

u = [u(N*dI-1:N) ;u(1:N+dI)]; x = x+dI*h; 
elseif dI >= 1 

u = [u(1:dI);u(dI+i:N)]; x = x+dI*h; 


end 


Moving the vanishing points to the side of the domain towards which the peak is travelling 
ensures that the periodicity of the solution is not destroyed. However, since we also require 
that any non-periodic solution, such as a soliton, be zero in the vicinity of the boundaries 
it is also possible to use a block of zeros for the new points without altering the solution. 
If the solution is not actually zero near the boundaries then this has the effect of filtering 
out the points on the vanishing edge of the computational domain. We will investigate any 
effects that using a moving computational domain may have on the accuracy of the solution 
in chapter []] once we have introduced higher order methods which we would like to use in 


conjunction with such a moving domain. 


An additional technique to deal with the periodic boundary conditions is to use absorbing 
boundaries to help ensure that outgoing waves which might be converted to incoming waves 


at the other side of the domain have little effect on the solution. 


3.2 Absorbing boundaries 


Because we are imposing a finite, implicitly periodic computational domain on a problem 
defined on the entire real line, it is necessary to have boundaries which absorb outgoing 
waves. For the NLSE this can be achieved by adding a complex potential term which is 
non-zero only in the vicinity of the boundaries to the right-hand side of equation (2.13) 
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Knapp, 1995|. We consider two possible potentials: 





vicos (3 TTL ) — vi, Хр<хж< 1, 








2Xr-aL 
Va = 4 0, ZL<EZ<ZER (3.2) 
vi COS (425) —Uj, TRIT XR 
and 
0; COS (rg) - vs Х„<ххж< TL 
Y= 4 0, ZL<T<ER (3.3) 
0; COS (т-а) -vi XREzXXg 


where X, and X р are the left and right boundaries of the computational domain, x; and £R 
are the left and right ends of the part of the computational domain where the computation is 
accurate and there is no absorption and v; is the (complex) amplitude of the potential. The 
shapes of the potentials were chosen to present a smooth profile to the outgoing waves so as 
to help reduce reflections. The magnitude of the potential v; can be tuned so as to minimise 
the sum of the amplitudes with as small a length of the absorbing region as possible. The 
width of the absorbing region, and the magnitude of the potential were selected by allowing 
a wave to “collide” with the boundary and then measuring the amplitudes of the reflected 
and transmitted components of the wave. Figures [3.2] and [8.3] show two such collisions for 
the two potentials with v; — 8i and using 2096 of the domain for the absorbing region. 
Also shown are the profiles along x = X, and along t = Т. These show the profiles of the 


transmitted and reflected waves respectively. 


'The results of measuring the magnitudes of the reflected and transmitted waves from 
a soliton of amplitude 1, for different choices of v; and using 2096 of the domain for the 
absorbing region are given in Table [3.1] 

We see from Table [3.I]that the two potentials treat the transmitted and reflected com- 
ponents of the wave differently. Although the potential V, is excellent at absorbing the 
transmitted component of the wave, it does not perform so well at preventing reflected 
components. Table B.I] also includes data for the numerical approximation of the integral 
of the magnitude of the solution at t = T' = 10 over the spatial domain calculated using 
Simpson's rule. This is a similar quantity to the remaining “mass” of the solution after the 
collision. This gives a second measure of the effectiveness of the absorbing boundaries and 
takes account of the fact that the profile of the wave may change after the collision meaning 
that the amplitude may not always be the best measure of the magnitude of the transmitted 
and reflected waves. Although the lowest total amplitude after the collision is given by the 


potential V, with v; = 15i, the measure given by the integral over |U(x,t = T)| is most 
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Figure 3.2: Collision of a soliton with the boundary of the computational domain with absorbing 
potentials, V, (top) and V, (bottom); v; = 8i, N = 128, q = 2. 
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Figure 3.3: Profiles of a soliton after a collision with the edge of the computational domain with 
the absorbing potentials, V, (left) and V, (right); v; = 81, N = 128, q = 2. The top figures show the 
profile along x = Xr; the lower figures, along t = T = 10. 
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2i | 2.9x 10-1 | 0.067 | 2.4 | 6.2 x 10-2 | 0.014 | 0.5 | 3.5 x 1071 | 7.6 x 107 
Ai | 1.0 х 107! | 0.013 | 0.7 | 8.3 x 107? | 0.022 | 0.7 | 1.2 x 107! | 3.0 x 107? 
6i | 4.6 x 107? | 0.013 | 0.4 | 1.7 x 107? | 0.030 | 1.0 | 5.9 x 1077 | 3.2 x 107? 
8i | 2.2 x 107? | 0.016 | 0.4 | 4.0 x 107^ | 0.036 | 1.3 | 3.7 x 1072 | 3.6 x 107? 
10i | 1.0 x 1072 | 0.018 | 0.5 | 1.1 x 107^ | 0.041 | 1.5 | 2.9 x 1077 | 4.1 x 107? 
12i | 5.6 x 107? | 0.020 | 0.6 | 3.7 x 107? | 0.045 | 1.8 | 2.6 x 1072 | 4.5 x 107? 
15i | 2.4 x 107? | 0.023 | 0.7 | 7.9 x 1076 | 0.053 | 2.0 | 2.5 x 1072 | 5.3 x 107? 
20i | 6.3 x 10-4 | 0.028 | 0.9 | 8.5 x 1079 | 0.062 | 2.4 | 2.9 x 107? | 6.2 x 107? 
25i | 1.8 x 10-4 | 0.032 | 1.1 | 8.7 x 10% | 0.070 | 2.8 | 3.3 x 107 | 7.0 x 107? 


Table 3.1: Measures of the effectiveness of the absorbing boundaries for V, and V, with different 
values of vj. вшуь = max, U(z = Xz,t)|, ешъ = шах, |U(z,t = T)| , Гь = J Ea |U(z, = T)|dz 
and Таль = Sa/b + Ea/b- 


favourable for the potential Va with v; between 6; and 8i. As a balance between these two 
measures we propose the potential V, with v; = 82 and using 10% of the computational 
domain on each boundary as a good choice to absorb any outgoing waves at the boundary; 
it is shown in Figure [3.4] 
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Figure 3.4: The complex potential for absorbing boundaries as given by (3.2) with X, = —25, 
Хр = 25, хр = —20, хв = 20 and vi = 8i. 


3.3 Mass conservation 


Since the soliton solutions of the NLSE are not periodic we no longer have the exact conser- 
vation of mass which we saw in the previous chapter for periodic solutions calculated with 
the Fourier method. The conservation properties of the numerical scheme are responsible for 
preserving much of the interesting behaviour which we expect from the analytic solution, for 
example, recurrence. Without conservation any perturbations of the numerical solution can 
grow without bound and one of the important analytic properties of the NLSE is lost. We 


therefore investigate the mass conservation property of the split-step Fourier methods with 
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a soliton solution and for the same initial conditions and parameters as used in Figure B.1] 
above but with the integration period extended to Т' = 40 in order to be more likely to see 
any effects which are slow to appear. Because of the longer integration period we make use of 
a moving spatial domain. The exact mass for these initial conditions is calculated according 
to equation (2.15) and gives M = 2. The mass of the numerical solution is calculated by 
approximating (2.15) with Simpson’s rule. The difference between the exact initial mass and 
that calculated from the discretised initial condition using Simpson’s rule is in the order of 
107°. 
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Figure 3.5: The percentage change in the total mass of the solution as a function of time. The 
solution was computed on a moving spatial domain with Т = 40, all other parameters are identical 
to those in Figure BI] 


Figure [3.5] shows that after a brief period of apparently no change in the total mass, the 
percentage error in the numerical value of the total mass begins to increase approximately 
linearly. The final percentage error in the mass was 2.57 x 107° indicating that, although 
mass is not conserved exactly, the change in mass is small enough that we can still claim 
to be approximating the true behaviour of the solution on such time scales. The percentage 
change in the mass corresponding to the calculation in Figure [3.T]is 1.07 x 1076, nearly four 


orders of magnitude smaller than the error of the same solution in the infinity norm. 
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3.4 Soliton collisions 


One of the interesting properties of soliton solutions is that they can exist in close proxim- 
ity and even overlap or collide without losing their profile, except momentarily during the 
collision. We see this if we begin with an initial condition given by two versions of equation 
(3.1) shifted in space and evaluated at t = 0. We take a linear superposition of two such 
solutions with ту = —8, v! = 0.3, ф = 0 and С! = 0.5 for the first and 22 = 7, v? = 0.6, 
62 = 0 and С? = —0.2 for the second. The parameters used for the integration are N = 512, 
q= 2, т = 0.01, L = 80 and T = 15. 

Figure B.ö]shows the numerical results for the collision. It is clear from both the plot of 
the magnitude of the solution and from the plot of the initial and final profiles of the solitons 
that the profiles of the solitons are maintained after the collision. The plot of the error in 
the mass conservation shows that as the two solitons begin to overlap the mass begins to 
fluctuate. These fluctuations grow as the collision proceeds but still remain very small overall 
(107796). Interactions between two or more solitons can present a computational challenge 
as the method of moving spatial domains is no longer applicable. However, as the solitons 
emerge from collisions with their profiles intact it would be possible to treat each separately 
on a moving domain when they are well separated and to test for impending collisions which 
would then be calculated on a single domain including both the solitons until they were 


again well separated. 
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3.4. Soliton collisions 
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Figure 3.6: A collision between two solitons, each with the general form given by (B.I) and with 
parameters zj = —8, v! = 0.3, ф = 0, (1 = 0.5, 22 = 7, 1? = 0.6, ф = 0 and (? = —0.2. The 
integration parameters are N = 512, q = 2, т = 0.01, L = 80 and T = 15. Below are shown the 
percentage change in total mass of the solution (left) and the profiles of the two solitons (right). The 
solid lines correspond to the profiles before the collision, the broken lines, afterwards. 





4: Higher Order Split-Step 
Methods 


4.1 Symmetric integrators 


Here we use the BCH formula to find symmetric splittings of higher order, first with exact 


coefficients, then later with coefficients given by a system of algebraic equations. 


4.1.1 Strang splitting and the symmetric BCH Formula 


In section [2.1] we used the BCH formula to find a unitary, first order approximation to 
the formal solution given by equation (2.3). We would like to be able to find higher order 
approximations to (2.3) which we can calculate numerically and which are still unitary and, 
therefore, still conserve the total mass of solutions of the NLSE. More generally speaking we 
wish to determine real coefficients су and d; such that for arbitrary, non-commuting matrices 


Aand B 
k 


exp(A(A + B)) = П exp(c;AA) exp(d;AB) + O(AP*1). (4.1) 
ј=1 

This gives an integrator of order р. For ће NLSE А = іт and, therefore, if the product in 
(LI) is symmetric in the coefficients c; and dj, then the integrator is clearly unitary. We have 
seen already that for the case p = 1 the BCH formula immediately gives exp(A(A + B)) = 
exp(AA) exp(AB) + О(А?), that is, су = dı = 1 and k = 1. The case for p = 2 is given by 
the well known Strang splitting (Strang, 1968}, also known as symmetric Trotter splitting 
[Trotter, 1959], with с = со = §, dj = 1, 4 = 0 and k = 2, thus 


exp(ACA Ву) e erp (5м) арар (5м) + O08). (4.2) 


Error bounds for the Strang splitting, with an unbounded operator A, are derived in 
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where it is shown that when applied to the NLSE with Fourier 
spatial discretisation, the error bounds for the Strang splitting depend on the temporal step 
size т and the smoothness of the initial condition и(2, 0), but not on the number of spatial 
grid points or on the number of time steps. 

In order to construct higher order symmetric splitting methods it is convenient to have a 
symmetric version of the BCH formula. We begin with the original form of the BCH formula 
as given by (1.44) which we multiply from the right by exp(AA) to get 


exp(AA) exp(AB) exp(AA) = exp(C(A)) exp(AA). 


We then rewrite C(A) as 





C(A) = AČ AI HE AT + 


where the C; are the coefficients of the BCH formula as given in (43). We are now able to 
apply the BCH formula (1.44); 

e484 = exp(AC) exp(AA) 

= exp (ле + А) + 


A 
+21 


3 
12 


A? 


5-8. AL 75 (Сб, А] + 14 A4, C1) 


[C, A, A,C] + ot) 





= exp ( € +20 t A?C3 + АСА + A) 


+2 ([C1, A] + [ACs, A] + [A?C3, A]) 





"12 ([C1, Ci, A]+[Ci, АС, АЈ АС, C, Al+[A, A, C1]+[A, A, AC3]) 


M 
+21161, 4,4,1] + oo) 


Setting in the expressions for C; i — 1,2,... and rearranging so that all the commutators 
can be written in the form [W, X,..., Y, Z] we find that the terms with even powers of А 
vanish to give 


1 


е^Ае^Ве^А = exp (лел BIER (5и, A, B] - c[B. B, Al) + oo) е (4.3) 


If we replace the argument A with A/2 so that both A and В are propagated over the 


same interval we get the familiar second order symmetric integrator (42) which we can write 
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as 





1 1 
p2(X) = exp (5м) exp(AB) exp (5м) = exp(AD, + à? D3 + °Ds ФАР; +--+) (4.4) 


where 
Dı = A+B 
Ds = d[B,B, A] - AIA 1A, В] 
Ds = адаан тунан 


+5014, В,В,В,А| + 355[B, 4, А, A, B] (4.5) 


sts 
– 15А, A, B, B, A] + |В, B, A, A, B] 


In order to see why (43) contains no even powers of A we observe that (A) is time- 
reversible with the inverse of exp (3AA) exp(AB) exp (3AA) obtained by changing the sign 
of A. The same must therefore be true for the right-hand side of (43) and EA) It follows 
that Dm = 0 for m = 2n, n € N. More generally, for any time-reversible integrator of the 
form (4.1) with the expansion 





P(A) = exp(Ay + A792 + Аз + Ауд s), 








we have 72 = y4 = у =+- = 0 |Yoshida, 1990]. To see this we form the product of (Л) 
and its inverse 
PX) = exp(-Ayı + Xa — Аз M — ++): 


From the BCH formula (LA) we see that e(A)g(—A) = exp(2A?72+O0(A3)). Since y(A)y(—A) 
= I, the argument of the exponential function, 2A?55 + O(A?), must vanish for all A hence, 
y2 = 0. The same product now gives y(A)y(—A) = exp(2A4y4 + О(А?)) which requires that 
Ya = 0. Continuing in this fashion we get yg = yg = ··· = 0. Therefore, any symmetric 


time-reversible integrator of the form (£I) has even order accuracy. 


4.1.2 Symmetric integrators with exact coefficients 


We are now able to use the symmetric form of the BCH formula to find exact values for 
the coefficients in (LI) so as to obtain higher order symmetric integrators. A fourth order 
integrator can be obtained from a symmetric repetition of the second order integrator (ДД) 


with the form 
Фа(А) = фә(ж1А)фә(жоА)ә (ж А), (4.6) 
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where the real unknowns х0 and x; must be determined such that the method has order four. 
Applying the symmetric BCH formula given by (43) to (46) with (21А) = exp(z1AD1 + 
23А D1 + хўА° D1 +--+) and фә(жоА) = exp(zoADı + 2А Dı + 3529 Di +--+) gives 


Фа(А) = exp(A(2z1 + zo) Di + A (221 + 23) Юз + ^°(2х? + 32) Ds +---). (4.7) 
In order for (7) to be a fourth order integrator we need 
27;--29 —1, and 22 +28 = 0 (4.8) 


so that (£7) becomes 
ya(A) = exp(A(A + B) + О(А?)). 


The unique real solution to (48) is 


23 1 


20 = 1 = ——- 
2— 23 


EC EE 
2— 23 


The integrator (4.7) can be written in the form (LT) with the coefficients di = d3 = 21, 
d2 = zo, сү = сд = 4/2 and сә = c3 = (10 + 21). This is the same as the fourth order 


integrator given in |Forest and Ruth, 1990}. 


Starting with the fourth order integrator given by (47) and continuing in the same 


fashion as above it is easy to obtain a sixth order integrator by using the symmetric product 
y6(A) = фла(у1А)л (уоА) а (из А). (4.9) 
The fourth order integrator has an expansion of the form 
ya(A) = exp( AD, + A Ds + A! D; + O(9)) (4.10) 
and so once again, applying the symmetric BCH formula (£3) to (£9) gives 
y6(A) = exp (Al2yı + yo) D + A? (2y? + 90) Ds + O(A7)) (4.11) 
which has order six when 2y; + yo =1 and 2y + Và — 0, or 


1 
25 1 


Yy = —— 15 yı = т. 
2 — 25 2— 25 


More generally, if we already have a symmetric integrator of order р = 2n, n = 1,2,... 
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given by фәһ (А) then we can determine the exact coefficients for a symmetric integrator of 


order p = 2n + 2. This is given by the product 


Q2n42(4) = ves (21А) Pan (202) 920 (21А). 
zı and zo must satisfy 221 + zo = 1 and gren + ze — 0 from which it follows that 


gm 1 
2 = — a5 — Hs 
2— 27 2— 27H 
We are therefore able to construct symmetric integrators of arbitrary (even) order from 


symmetric compositions of the second order symmetric integrator pa(A). 


By rewriting the higher order order integrators in terms of the second order integrator we 
can determine the values of the coefficients c; and d; in { Т). For the sixth order integrator 
(ZIT) given above we have 

Фв(А) = YalzıyıA)p2(zoyıA)Pp2(zıyıA) х PalzıyoA)p2(zoyoA)P2(XıYyoA) 
x фо(ату А) о (хо А) о (2101 А) 
which gives dı = da = d; = dg = zıyı, da = dg = zoyı, da = dg = x1yo, ds = х0уо, and 
сү = Fd, Cig = йә, 61 id), i= 2,3,...,9. 








Although it is possible to develop symmetric methods of arbitrarily high order with 
relative ease in this fashion, such methods require the p2(A) integrator 3"~! times to attain 
order 2n. This means that the number of exponential operators, k = 3", grows rapidly as 
the order of the integrator increases. We next look for symmetric integrators which require 


fewer exponential operators. 


4.1.3 Symmetric integrators with fewer steps 


Rather than constructing an integrator of order 2n + 2 recursively from repetitions of an 
order 2n integrator, it is possible to instead begin directly with a symmetric product formed 
from the second order integrator (4.4). Here we give an outline of the method of Yoshida 
which allows us to find symmetric integrators which use fewer y2(A) opera- 
tors. We define an integrator formed from a symmetric product of m — 2n -- 1 second order 


operators 


v" (A) = ga(wnA) х +++ x qa(wiA) x qga(woA) x PalwıA) x +++ x Yo(WndA), (4.12) 
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where wo,..., Wn are real coefficients to be determined. Using the symmetric BCH formula 


(£3) we are able to express (4.12) as a single exponential function similar to (27); 


Ф"(А) = exp(Ad(1,n)D, + Аа(3, п) Юз + A (d(5,n)Ds + e(5,n) Es) 
+X (d(7, n)Dr + e(7, n)Er + FCT, n)Fr ЕЯ 9(7, n)G7) + О(А)), 
(4.13) 


where we have substituted E5 = [2у, Dj, D3], Е; = |D1, D1, D5], Е; = [D3, D3, Dı] and 
G7 = [D4, Dj, Dı, D, Рз], and where the d(i, j), e(i,j), (0,3), and g(i,j) are nonlinear 
combinations of the unknowns wo,...,Wn. By applying the symmetric BCH formula (4.3) 


to 
Ф?Ү!(Х) = ea(us i A) p Qa (uma) 


and comparing both sides one obtains a recursion for the coefficients d(i, j), eli, 7), f(i, J), 
and g(i, j), along with the initial values, in terms of wo. This gives a nonlinear simultaneous 
system of equations which must be solved for wo,...,Wn. In order that the integrator uses 
the minimal number of exponential operators and has order six п = З is necessary and 
sufficient, while for an eighth order integrator we need n = 7. It is always possible to 
reduce the order of the system of equations by one through setting a(1,n) — 1 which gives 


wo = 1— 2(0 +w2 +: Бо). The resulting set of algebraic equations has the general form 


ti (Was, .-, Wn) = 0, 
yo(w1, W2,...,Wn) = 0, 
Yn(W1, We,...,Wr) = 0. 


and must be solved numerically. 


Tables (ZI) and (42) list some of the possible solutions given in [Yoshida, 1990] for 


orders six and eight respectively, and where the system of algebraic equations was solved 
using the Brent method |Brent, 1973}|Brent, 2002|. This method has the advantage that it 
does not require the Jacobian matrix Oy;/Ow; which, for the function above, is difficult to 


derive. 


These methods require 21 exponential functions for order six and 45 for order eight. In 
contrast, the methods from section A.1.2]require 27 exponentials at order six and 81 at order 
eight, hence, the methods given by (4.12) with the coefficients from Tables ÆI] and [4.3] offer 


4.1. Symmetric integrators 71 


Solution A Solution B 
wj —0.117 767 998 417 887 x 10! —0.213 228 522 200 144 x 10! 


wə 0.235 573 213 359 357 х 10° 0.426 068 187 079 180 х 1072 
w3 0.784 513 610 477 560 x 10° 0.143 984 816 797 687 x 10! 


Solution C 
wi —0.214 403 531 630 539 x 10! 


we 0.152 886 228 424 992 x 1072 
шз 0.144 778 256 239 930 x 10! 


Table 4.1: Yoshida’s coefficients for possible sixth order symmetric integrators of the form { Т2). 


significant computational savings. 


In the cases above we chose m = 2n + 1 such that (4.12) used the least number of фә 
operators possible for the desired order. It is however possible to use larger values of m 
which leads to free parameters in the system of algebraic equations. The free parameters 
can then be chosen to ‘tune’ the integrator so as to give the smallest possible coefficient 
for the leading error term. For integrators of the form (£12) m is always odd, however, by 
omitting the central yo(woA) stage from (£12) it is also possible to obtain methods with 
even m. McLachlan [McLachlan, 1995], however, found no cases where it was advantageous 


to take m even as doing so only provides the same number of unknowns as for m — 1. 


For order four with m = 5 one can improve on the m = 3 method of section [I.1.2] Suzuki, 
Suzuki, 1990], gives an exact coefficient method with шу = wa = 1/(4 — W4). McLachlan 


McLachlan, 1995] improves on this slightly; the necessary coefficients are given in Table 3] 


Numerical searches indicate that for order six with m = 7 Yoshida's coefficients given 
as Solution A in Table ET] yield the optimal error coefficient. For m = 9 and m = 11 
methods with lower leading error coefficients are available. Both m = 9 and m = 11 yield 
approximately the same error coefficient and so the method with m = 9 as given in Table 
[4.3] is to be preferred (McLachlan, 1995]. For order eight and m = 15 the method with the 
coefficients given by Solution D in Table [4.2] is the best of the the five possibilities given 
by Yoshida, however, a numerical search by McLachlan found 100 possible solutions with 
‘reasonable’ parameters, the best of which has a leading error coefficient about an order 
of magnitude smaller than that of Yoshida’s method. The same result is also given (to 25 
significant figures) in Suzuki, 1994]. Taking m = 17, an error coefficient another order of 


magnitude lower was found. For m = 19 the best error coefficient was marginally larger 


McLachlan, 1995]. 
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Solution A Solution B 
wı  —0.161 582 374 150 097 x 10° —0.169 248 587 770 116 x 1072 
t9  —0.244 699 182 370 524 x 10° 0.289 195 744 315 849 x 10° 
из —0.716 989 419 708 120 x 107? 0.378 039 588 360 192 x 1072 
w4 0.244 002 732 616 735 x 10+ —0.289 688 250 328 827 x 10" 
Ws 0.157 739 928 123 617 x 10° 0.289 105 148 970 595 x 10" 
We 0.182 020 630 970 714 x 10° —0.233 864 815 101 033 x 10° 
Шт 0.104 242 620 869 991 х 10° 0.148 819 229 202 922 x 10° 
Solution C Solution D 
w 0.311 790 812 418 427 x 10° 0.102 799 849 391 985 x 10° 


wo —0.155 946 803 821 447 x 10! —0.196 061 023 297 549 x 10° 
ws —0.167 896 928 259 640 x 10! 0.193 813 913 762 276 x 10° 
wa 0.166 335 809 963 315 x 10! —0.158 240 635 368 243 х 10° 
ws —0.106 458 714 789 183 x 10! —0.144 485 223 686 048 x 10° 
we 0.136 934 946 416 871 x 10! 0.253 693 336 566 229 x 10° 
w7 0.629 030 650 210 433 x 10° 0.914 844 246 229 740 x 10° 





Solution E 
WI 0.227 738 840 094 906 x 10-1 


wa 0.252 778 927 322 839 x 10! 
w3 —0.719 180 053 552 772 x 1071 
wa 0.536 018 921 307 285 x 107? 
ws  —0.204 809 795 887 393 x 10! 
wg 0.107 990 467 703 699 x 10° 
ит 0.130 300 165 760 014 x 10! 


Table 4.2: Yoshida's coefficients for possible eighth order symmetric integrators of the form (£12). 


4.2 Numerical results 


Here we investigate the order behaviour and conservation properties of the split-step Fourier 
method over a range of orders and when applied to soliton solutions. We mention again 
that the Fourier method assumes a periodic solution. It is therefore important that we 
consider the effect of the spatial discretisation on the theoretical order behaviour and mass 
conservation. Here we consider methods using first order (Lie) and second order (Strang) 
splitting, as well as two methods each of orders four, six and eight. The coefficients of the 
methods are given in Table [4.3] the order four methods are the exact coefficient methods 


with three or five stages. The methods of order six use seven or nine stages and the order 
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4 Up. 
wz = 1/(4- V4) 
0.625 466 428 467 670 01 шә = 0.28 
6 —1.177 679 984 178 871 006 95 ш» = 0.235 573 213 359 358 133 68 
0.784 513 610 477 557 263 82 
—0.843 265 623 387 734 608 55 ш» = 0.129 466 948 913 475 358 06 
0.555 497 023 712 478 39916 w4 = 0.1867 
8 = 0.315 293 092 396 766 596 63 ws = 0.334 624 918 245 298 183 78 


= 0.299 064 181 303 655 923 84 w4 =—0.573 862 471 116 082 266 66 
= 0.190 754 710 296 238 379 95 we =—0.409 100 825 800 031 594 00 
0.741 670 364 350 612 953 45 

= 0.314 956 683 916 294 857 89 wo = 0.278 335 500 393 679 651 31 
= 0.144 405 941 080 012 041 06 w4 =—0.409 552 343 420 851 419 34 
= 0.185 146 935 716 587 732 65 wg ——0.410 175 371 469 850 137 53 
= 0.581 514 087 105 250 962 43 ws = 25/194 





Table 4.3: Coefficients for symmetric methods composed of symmetric products of фә with optimal 
leading error coefficients. The coefficients for wo can be determined by wo = 1 — BEN w;. All 
coefficients are exact to 20 significant figures. 


eight methods use 15 or 17. All calculations were run on a moving computational domain 
of width L = 50, both with and without the absorbing boundaries described in section [3.2] 
In order to see any effects of evaluating the nonlinear term in the middle of the symmetric 
step, as opposed to evaluating it at the beginning and end of the symmetric step, we present 
results for both cases. We denote the two cases by LNL and NLN respectively! 

The problem solved is the same as the single soliton problem in chapter B] with initial 


conditions given by 
u(z,0) = exp (i (22 + 5)) sech(x). (4.14) 


'The exact solution, 
u(x,t) = exp (i (2 — 3t + =)) sech(x — At), (4.15) 


is used as a comparison. Using д = 2 we integrated the NLSE to Т = 2 for a range of 
temporal and spatial step sizes. The total error at the end of the integration interval was 


calculated in the La norm, according to 


1 
2 


М 
err = |h S "|U(zj T) -u(zj, T^ |, (4.16) 
j=l 


! We will sometimes neglect to write LNL, i.e. if a method is not specified as being NLN then it is LNL. 
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where h is the spatial step size and U(r;, Т) and u(x;,T) are the numerical and exact 
solutions at x; = jh respectively. 
4.2.1 Accuracy of Lie and Strang methods 


Figure[A.T]shows the error plotted as a function of the temporal step size for the Lie splitting 
methods with h = 50/128, 50/256, 50/512 and 50/1024. 
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Figure 4.1: Error as a function of step size for the Lie splitting scheme for 128, 256, 512 and 1024 
spatial grid points. The solid and the dotted lines correspond to results calculated with and without 
absorbing boundaries (ab) respectively. The thick solid line shows a line of gradient one, i.e. the 
theoretical order behaviour of the first order scheme. 


From this we see that the error follows the expected behaviour for a first order method. 
The error of the method is essentially unaffected by the number of Fourier nodes used in the 
spatial discretisation with the exception of the results for N = 128 when small (т < 1073) 
step sizes are used. For this case we see that the decrease in the error obtained by decreasing 
the integration step size begins to slow slightly once the error reaches about 7 x 1074. We 
suggest (and later show) that this effect is due to the spatial discretisation error which gives 
a lower bound for the minimum error of a method irregardless of the integration step size. 


For the Lie splitting methods the absorbing boundaries do not seem to have a noticeable 
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effect on the final error. 
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Figure 4.2: Error as a function of step size for the Strang (LNL) splitting scheme for 128, 256, 512 
and 1024 spatial grid points. The solid and the dotted lines correspond to results calculated with 
and without absorbing boundaries (ab) respectively. The thick solid line shows a line of gradient two. 
The results for the NLN scheme are identical. 


The results of the Strang splitting method (LNL) shown in Figure [4.2] are essentially 
similar to those of the Lie splitting with the error following the expected behaviour of a 
second order method. Here we can see that for N — 128 the error remains constant at about 
5 x 107^ when the temporal step size is decreased below т = 1/100. In this case the error 
curve plateaus off indicating that a spatial step size of h = 50/128 is too coarse and results 
in the spatial discretisation error limiting the total error. The results for the Strang splitting 
method with NLN are essentially, (i.e. to graphical precision), identical to those of the LNL 
method and hence we omit them. Before we consider the accuracy of higher order methods 


we investigate the effect of spatial discretisation in more detail. 


4.2.2 Effect of spatial discretisation 


The extent of the spatial discretisation error can be quantified by comparing the numerical 


approximation of the second derivative with its exact value. The second derivative of the 
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solution (£15) at t = 0 is given by 


9? 
a: ue t) = — exp(i(2x + m/2))sech(x) (4 + 4i tanh? (x) + sech”(z)) А 

е t=0 
We can compare this with the numerical approximation to the second derivative by using 
the DFT to calculate the Fourier coefficients йк, k = —N/2,..., N/1 — 1 of the discrete 
solution of length N given by the discrete version of the initial condition (£14). We then 
use the second derivative operator in Fourier space along with the trigonometric interpolation 


polynomial 


ща) рз) = У) dette 


to give а truncated Fourier series representation for the second derivative of the initial con- 


dition. This is given by the trigonometric interpolating polynomial for the second derivative: 





д? dp X mb. irem 
az ula t) cB = ‘> (-4 72 а) gene (4.17) 


t=0 р 


То measure the error іп the spatial discretisation of the second derivative we calculate 
the L5 norm of the difference between the exact second derivative and its numerical approx- 
imation calculated via equation (4.17) and evaluated on a fine mesh. We scale the norm by 
the square root of the number of grid points used on the fine mesh. The results calculated 
on the interval [—25,25] for 128, 256, 512 and 1024 Fourier nodes are given in Table 4 
Different values for m, the number of grid points on the fine mesh, are used to confirm that 
the mesh is fine enough to give a true indication of the error over the continuous interval 
[—25, 25]. Rather than evaluating the sum in equation { ТД) directly for each x on the fine 
mesh we use the FFT and IFFT to interpolate the values for the second derivative once 
they have been calculated using the Fourier method in the usual fashion. We see from these 
results that for N = 128 Fourier nodes the error in the numerical approximation to the 
second derivative is approximately 107. This gives an indication of the size of the error in 
the spatial discretisation at each step of the splitting method and confirms that the lower 
bound of approximately 5 х 1074 in the total error for the Strang splitting with N = 128 
spatial grid points is most likely to be due to the spatial discretisation error. The results 
also predict that we will see similar behaviour for N = 256, 512 or 1024 Fourier nodes once 
the total error is in the vicinity of 1078 to 10—19. Because attaining such an accuracy with 


a method of order one or two would require very many more time steps and a much smaller 
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m Lo error N m Lo error N 

2048 | 1.53 x 107 8192 | 1.53 x 107 
1.99 x 1078 1.99 x 1078 
1.63 x 10—10 1.63 x 10-10 
4.59 x 10—10 4.58 x 10-10 

4096 | 1.53 x 107 128 16384 | 1.29 x 1079 | 2048 
1.99 x 1078 | 256 16384 | 3.66 x 107? | 4096 
1.63 x 10-10 | 512 16384 | 1.03 x 1078 | 8192 


4.58 x 10—10 | 1024 


Table 4.4: Error in the spatial discretisation of the second derivative evaluated on the interval 
[-25,25] with N = 128, 256, 512 and 1024 Fourier nodes. The number of points used for the 
fine mesh is indicated by m. 


integration step size introducing the possibility of rounding errors, we delay the investigation 
of this claim until we look at the accuracy of the higher order splitting methods in section 
[4.2.4] Although increasing the number of Fourier nodes from N = 256 to N = 512 results 
in a decrease of roughly two orders of magnitude in the error in the approximation to the 
second derivative, increasing the number of Fourier nodes further to N = 1024 does not 
give any increase in the accuracy of the approximation. Rather the error for N = 1024 is 
slightly larger than that for N — 512 in all cases. This suggests that rounding errors limit 
the accuracy of the approximation to about 10-1 and that increasing N above the smallest 
integer at which this accuracy can be obtained has the effect of increasing the error slightly 
due to the accumulation of these errors. This is supported by the results for N = 2048, 4096 
and 8192 Fourier nodes which still give errors of about 107° in the approximation to the 
second derivative. We therefore use a minimum of 256 spatial grid points and a maximum 


of 1024 in the calculations for the higher order methods. 


4.2.3 Effect of a moving computational domain 


In section [3.1] we introduced the idea of a moving computational domain which could be 
used to ensure that the peak of the solution remained at the centre of the spatial domain. 
In order to see if this has any effect on the accuracy of the final solution we first compute 
a solution using a moving domain so as to find the final values of the spatial domain. We 
then create a fixed domain which covers the whole spatial domain spanned by the moving 
domain during the integration period. We calculate the solution on this fixed domain using 
the same spatial and temporal step sizes as we did for the moving domain and with initial 


conditions generated by the same function (i.e., equation (4.14) in both cases). The exact 
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solution as given by equation (£15) can then be used to calculate the error in the solutions 
corresponding to the two different domains. 

For the case of a moving domain initially given by [—25,25] and using N = 256 (i.e., 
h = 50/256) we find that for an integration period of T' = 10 and a time step size of 
T — 0.02 the computational domain shifts by 204 points during the course of the integration. 
The corresponding fixed domain is therefore given by a domain of 460 points beginning at 
x = —25 and increasing with a spatial step size of h = 50/256. 


We use two different measures of the error of the solution for the fixed domain: 


1. We calculate the error by using equation (£16) to compare the numerical and the 


exact solutions at all points in the fixed domain. We will denote this error err fe 


2. We calculate the error by using equation (416) to compare the numerical and the 
exact solutions at only those points in the fixed domain which correspond to points in 


the moving domain at the end of the integration period. That is, we use only the last 
256 


fcc 


256 points of the 460 points in the fixed domain. We denote this error err 


The reason for using the two different measures of the error is that, although both calculations 
use the same spatial step size, the calculation on the fixed domain uses a greater total number 
of Fourier nodes due to the extra length of the computational domain. The initial condition 
is therefore better sampled and the resulting Fourier representation of the solution includes 
higher frequency terms than the corresponding solution calculated on the moving domain. If 
we therefore compare the errors in the solutions only using those points on the fixed domain 
which correspond to the moving domain we are effectively comparing a better sampled 
solution at only some of its points and and might therefore expect a more accurate result. 
On the other hand, if we use the error calculated on the entire fixed domain then we are 
comparing errors which correspond to different spatial domains sampled with the same grid 
spacings and which may not therefore be well matched. By using both these measures of 
the error in the fixed domain solution we hope to overcome, or at least unveil, any such 
problems. 

The error in the solution from the moving domain is calculated as usual using equation 
(4.16) and is denoted errmy. 

The errors were calculated on the two domains using a range of splitting methods from 
order two through to eight. The results are given in Table [1.5] 

The results in Table [.5]show that any effect due to the use of the moving domain is in 
the order of 1071}, at or below the level of the spatial discretisation error and well below 
the error in all but the most accurate results. The measure of the error given by ero is 


consistently less than (or equal to up to 4 significant figures) the error of the moving domain 


4.2. Numerical results 79 


method 8,17 8,15 6,9 6,7 





етте» | 1.273 x 10710 1.281 x 10710 6.277%10- 7.406 x 107° 
err? | 8.799 х 1071! 8.762 х 1071 6.245 x 10719 7.404 x 107? 
errs” | 1.013 x 107'° 1.009 x 10719 6.307 x 10719 — 7.477 x 107? 


method 4,5 4,3 Strang 





етт» | 4.357 x 1077 1.476 x 10-9 4.077 x 1073 
erro | 4.357 х 10-7 1.476 х 1075 4.077 x 107? 
етт} | 4.357 х 10-7 1.476 х 1075 4.077 x 10? 


Table 4.5: Comparison of errors on fixed and moving domains for Т = 10, т = 0.02 and h = 50/256. 
The moving domain is initially [-25,25] and uses 256 grid points. The error on the moving domain 
is given by erra, via eqn. (£16). The fixed domain uses 460 grid points and covers all the points 


spanned by the moving domain during the course of the integration. The error over all 460 points is 


given by are the error over the last 256 points by err?}°. The methods used for the integration 


are classified by the pair p,s where p is the order of the method and s is the number of symmetric 
stages used. 


еттт». This is most likely to be due to the reason mentioned previously; we have used a 
larger number of Fourier nodes to sample the initial condition and have then used only some 


of these to calculate the error. 


4.2.4 Accuracy of higher order methods 


Figures [4.3] to 1] show the order behaviour for the methods of orders four, six and eight. 
In general, the plots show that the methods display the expected time integration order 
behaviour with the slopes of the error versus step size curves roughly matching the gradi- 
ents of the curves for the theoretical order behaviour. For the methods with no absorbing 
boundaries we see that the size of the error at which the curves flatten off agree well with 
the predictions of the estimated spatial discretisation errors. For N = 256 Fourier nodes 
this is about 1 x 10—10 while for N = 512 and N = 1024 the lower bound observed for the 
total errors is approximately 7 х 1071 in both cases. 

When absorbing boundaries are used with the methods the lower bounds for the total 
errors are one to one and a half orders of magnitude higher. For the methods of orders six 
and eight we see a very slight improvement in the accuracy of the method when using the 
absorbing boundaries. This suggests that the computational domain used in the calculations 
was sufficiently wide that the periodic boundary conditions had a negligible effect on the 
solution. 


It is possible that the difference between the accuracy of the methods with and without 
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Figure 4.3: Error behaviour of splitting methods of order 4, 6 and 8 applied to the NLSE using 256 
spatial grid points. Methods are classified using the pair (p,m), where p is the order of the splitting 
and m is the number of symmetric stages used. Results are given with (solid lines) and without 
(dotted lines) absorbing boundaries (ab) and for both the NLN (upper plot) and LNL (lower plot) 
methods. The thick broken lines are of slope 4, 6 and 8. 
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Figure 4.4: Error behaviour of splitting methods of order 4, 6 and 8 applied to the NLSE using 512 
spatial grid points. Methods are classified using the pair (p,m), where p is the order of the splitting 
and m is the number of symmetric stages used. Results are given with (solid lines) and without 
(dotted lines) absorbing boundaries (ab) and for both the NLN (upper plot) and LNL (lower plot) 
methods. The thick broken lines are of slope 4, 6 and 8. 
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Figure 4.5: Error behaviour of splitting methods of order 4, 6 and 8 applied to the NLSE using 256 
spatial grid points. Methods are classified using the pair (p,m), where p is the order of the splitting 
and m is the number of symmetric stages used. Results are given with (solid lines) and without 
(dotted lines) absorbing boundaries (ab) and for both the NLN (upper plot) and LNL (lower plot) 
methods. The thick broken lines are of slope 4, 6 and 8. 
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the absorbing boundaries would be more noticeable if the width of the computational was 
badly chosen. However, it seems that by using a computational domain such that the 
magnitude of the solution is near zero in the vicinity of the boundaries the need for absorbing 
boundaries can be avoided and a lower level for the bound on the maximum accuracy of the 
method can be obtained. For the fourth order methods no benefit from using the absorbing 
boundaries was apparent. 

The plots of the errors for the higher order methods also show the benefits which can be 
obtained by using a larger number of symmetric stages for a method of the same order so as 
to ‘tune’ the method and reduce the leading error term. For the fourth order methods we 
see a gain in accuracy of about one and a half orders of magnitude. The benefit decreases 
with increasing order of the methods; for the eighth order methods slightly under half an 


order of magnitude can be gained by using two extra symmetric stages. 


4.2.5 Mass conservation 


Using a Fourier method for the spatial discretisation means that the mass conservation of 
the splitting scheme is exact only for solutions which are periodic in space. Therefore, we 
now study the mass conservation properties of the schemes considered above for the case 
when the solution is a soliton and hence non-periodic. The error in the mass conservation 
is calculated by comparing the mass of the initial solution with the mass at time Т. As in 
chapter B] this is calculated by using Simpson's rule to approximate the integral of |U|? over 
the spatial domain. 

Figureff.ö]shows the error in the mass conservation for the Lie splitting scheme. From this 
we see that, although the mass is not conserved exactly, the error in the mass conservation is 
far smaller than the Lə error in the solution as given by (4.16) and shown in Figure LI] For 
T — 0.1 the error in the mass conservation (when absorbing boundaries are not used) is in the 
order of 1078 while the total error is close to 1071. Additionally, we see that although the 
Lie splitting is first order accurate, the error in the mass decreases with decreasing temporal 
step size at twice this rate, i.e. like a second order method. As in Figure ÆI] we also see 
flattening off of the curves corresponding to N = 128 due to the spatial discretisation error. 
We note that although effect of the absorbing boundaries on the solution was well below the 
error and hence not observable for the Lie splitting, the error in the mass conservation for the 
same method is small enough that the effect of the absorbing boundaries becomes apparent 
with the error in the mass conservation being approximately one order of magnitude greater 
when the absorbing boundaries are used. 

Similar behaviour can be seen in Figure [4.7] where the mass conservation properties of 


the Strang splitting methods are shown. 
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Figure 4.6: Error in the mass conservation of the Lie splitting scheme for the NLSE with 128, 256, 
512 and 1024 spatial grid points. The solid and the dotted lines correspond to results calculated with 
and without absorbing boundaries (ab) respectively. The thick solid lines are of slope one and two. 


Again, we see that the absorbing boundaries increase the error in the mass conservation 
by roughly an order of magnitude and that the error in the mass conservation decreases 
with decreasing step size at roughly twice the rate of the total error. This phenomenon, 
combined with the small size of mass conservation error even for larger step sizes, means 
that the error in the mass conservation quickly approaches machine precision. Before this 
becomes an effect, however, the limit for the accuracy of the spatial discretisation prevents 
the energy error from decreasing further. We see this occur when the mass conservation 
error is in the vicinity of 10714 and the integration step size is about 5 x 107%. Below this 
step size the accumulation of errors, probably due to both the spatial discretisation and to 
rounding errors, causes the error in the mass conservation to slowly increase again. 

Figures ES] to E. 10] show the mass conservation results for the higher order splitting 
schemes. As for the Lie and Strang splittings in Figures [4.6] and 7 we see that the error 
in the mass conservation is significantly smaller than the La error given by equation (4.16) 


and shown in Figures [4.3] to 1.5] 


The mass conservation error is lowest when the integration step size is such that the 
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Figure 4.7: Error in the mass conservation of the Strang splitting scheme with NLN (top) and LNL 
(bottom) ordering and both with and without absorbing boundaries (ab) for 128, 256, 512 and 1024 
spatial grid points. The thick solid lines are of gradient two and four. 
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Figure 4.8: Error in mass conservation for splitting methods of orders 4, 6 and 8 applied to the 
NLSE with (solid lines) and without (dotted lines) absorbing boundaries (ab). The spatial step size 
is h = 50/256. 


4.2. Numerical results 87 


























512, N-L-N 
107 
©: 4,3 
o 4,5 
-8 x: 6,7 
1 = , 
g ж: 6,9 
А. 8,15 
9 9 8,7 
10 9 -© 43- ab 
—- 4,5 – ab 
—- 6,7 – аб 
10719 H + 6,9 – ab 
—- 8,15 – ab 
геа =- 8,17 -ab 
e 10" L 
[0] 
d 
d 10°" 
E 
to? 
то. 
1079 = 
107'® А 
10 
step size 
512, L-N-L 
107 
©: 4,3 
0: 4,5 
-8 x: 6,7 
10 E.» 69 
А. 8,15 
Vv 8,17 
-9 || -© 43- ab 
10 E -= 45 ab 
эе 6,7 - ab 
— 6,9 - ab 
1079. + 8,15 – ab 
=- 8,17 -ab 











mass error 
© 











-4 


step size 


Figure 4.9: Error in mass conservation for splitting methods of orders 4, 6 and 8 applied to the 
NLSE with (solid lines) and without (dotted lines) absorbing boundaries (ab). The spatial step size 
is h = 50/512. 
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Figure 4.10: Error in mass conservation for splitting methods of orders 4, 6 and 8 applied to the 
NLSE with (solid lines) and without (dotted lines) absorbing boundaries (ab). The spatial step size 


is h — 50/1024. 
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total error in the solution has just reached the lower bound in the error due to the spatial 
discretisation. As was the case for the Strang splitting, this occurs when the mass conserva- 
tion error is in the range 10714 to 10713. For step sizes significantly smaller than this we see 
that the accumulation of the spatial discretisation (and possibly also rounding) errors has 
caused the mass conservation error to increase nearly to the magnitude of the total error. 
'This accumulation of errors shows the importance of ensuring that the total error in the 
solution remains above the lower bound due to the spatial discretisation error. 

We briefly note that the benefit of using extra symmetric stages to ‘tune’ the accuracy of 


the method is slightly greater for the mass conservation error than for the total error given 


by equation (4.16). 


4.2.6 Efficiency: work versus precision 


In order to make a judgement of the desirability of any particular method it is necessary to 
also consider the cost of the method relative to its performance. Figures EIT]to I. T4] show 
the work-precision diagrams corresponding to the above calculations. All calculations were 
performed on an Intel Pentium-4 2.8GHz PC with 1GB of RAM running MatLab 6.5. 

Figure[4.1]]shows a precision-work plot for the Lie splitting and the LNL Strang splitting 
results. For clarity, and because they are essentially identical to those of the LNL case, the 
NLN Strang splitting results are omitted. 

From Figure [TT] we see that for any reasonable error (i.e. less than about 1071) the 
second order Strang splitting methods are to be preferred over their first order counterparts. 
Although the efficiency of the methods is essentially unchanged by the use of the absorbing 
boundaries, we have seen in section [1.2.5] that such boundaries do decrease the accuracy of 
the mass conservation of the method. Since the mass conservation property of the NLSE 
is an important physical property which we would like to preserve exactly, the methods 
without the absorbing boundaries are to be preferred. For total errors larger than about 
107? using N — 128 Fourier nodes for the spatial domain is sufficient and results in the most 
efficient method. For total errors smaller than 1073, N = 256 is to be preferred. Increasing 
the number of spatial grid points above N — 256 gives no appreciable benefit. 

When the Strang splitting methods are compared with the higher order methods shown 
in Figures [I.12]to [4.T4]it is found that the fourth order methods are to be preferred over the 
Strang splitting methods. In particular, the order four method with five symmetric stages is 
optimal among the methods considered for errors down to about 1079. Beyond this point the 
preferred method is that of order six with nine symmetric stages. As for the case with the 
Lie and Strang splittings, the methods without absorbing boundaries are to be preferred. In 


addition to giving worse results for the conservation of mass and increasing the lower bound 
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Figure 4.11: Precision-work plot for the Lie (dotted lines) and the LNL Strang (solid lines) splittings. 
'The results for the NLN Strang method are essentially identical to those for the LNL Strang method. 


in the total error, the absorbing boundaries result in any marginal increase in accuracy being 


negated from the point of view of efficiency by the extra cost in applying such boundaries. 


We see that in general the methods with additional symmetric stages are more efficient 
than their counterparts with fewer symmetric stages but of the same order. The point at 
which the eighth order methods appear about to become more efficient than those of order 
six occurs at the same point at which the errors reach the lower bound due to the spatial 
discretisation; about 10-19. This makes the eighth order methods less efficient than the 


methods of order six, at least when calculations are performed in double precision. This 


agrees with the findings of Bandrauk and Shen |Bandrauk and Shen, 1994]. 
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Figures [4.12] to [4.14] also show that the effect of increasing the number of spatial grid 
points above N = 256 is only to shift the precision-work plots to the right, that is, to increase 
the cost of the computations with no appreciable increase in accuracy. This highlights the 
importance of correctly matching the number of spatial points to the domain/problem. We 
also note that the LNL methods are more efficient than their NLN counterparts; an indication 
of the fact that the nonlinear term is slightly more expensive to evaluate than the linear term. 
This is most pronounced for the case of N = 256. 

In this section we have presented results for a variety of different implementations of 
the split-step Fourier method applied to the NLSE for a range of orders and step sizes. 
In doing so we have illustrated the various properties of the different schemes. All the 
results obtained in this section have been compared with the known exact solution of the 
NLSE. An alternative technique would be to calculate a reference solution with a reliable, 
standard time integrator of known accuracy. Doing so, however, would remove the effect 
of the spatial discretisation from the comparison. By comparing the numerical results with 
exact solutions and considering the effects of spatial discretisation separately, we are able 
to study the interaction between the spatial and temporal errors. In this way we are able 
to suggest choices of methods and parameters which are likely to be well suited to related 
problems. An example of such a problem: by adding a potential term containing random 
steps to the right-hand side of the NLSE one is able to model the effect which joins and 
inhomogeneities in an optical fibre have on a pulse propagating along it [Knapp, 1995]. 

An additional reason for using the exact solution for comparison of results is that the 
calculations here were all performed on a shifting computational domain making an appro- 
priate comparison with a numerical reference solution from another method difficult as the 
shifts in the domain can depend on the accuracy of the method. 

We have seen here that the effect of the spatial discretisation error is to create a lower 
bound below which the total error in the solution can not decrease, irregardless of the 
integration step size. For the problem considered here this gave a bound of about 5 x 1074 
when using N = 128 and of about 10-1? for N = 256, 512 or 1024. Increasing the number 
of Fourier nodes was not able to reduce this bound. 

The absorbing boundaries introduced in chapter B] were found to have little effect on 
the accuracy of the solution and to have a negative effect on the mass conservation of the 
methods. It is therefore preferable, where possible, to avoid such artificial boundaries and 
instead use a computational domain which ensures that the solution is zero in the vicinity 
of the boundaries. 

For non-periodic solutions of the NLSE the splitting methods are no longer able to ensure 


that the mass of the solution is exactly conserved. It is possible, however, to ensure that the 


92 Chapter 4. Higher Order Split-Step Methods 

































































256, N-L-N 
10” T T 
© 4,3 
|. 4,5 
-3 x: 6,7 
10? b * 6,9 E 
A. 8,15 
у 8,17 
4 -© 4,3 – ар 
10 F = 45-ab ff 
— 6,7 - ab 
— 6,9 – ab 
10° L -A 8,15 - ab || 
-V 8,17 - ab 
m 
O is 
E 10° + 4 
ш 
107 $ Е 
10° $ E 
10° b XR 4 
10" 1 M ege esa Rec ge iul 
107 107 10° 10' 10 
Work — seconds 
256, L-N-L 
10° т І 
© 4,3 
a 45 
x. 6,7 d 
ж: 6,9 
А. 8,15 
V. 8,17 H 
-© 43-ab 
— 4,5 - ab 
—-6,7-ab H 
— 6,9 - ab 
-A 8,15 - ab 
-7- 8,17 -ab || 
п 
о 
pum c 
эм 
ш 
ASK 4 
A ; 
101° 1 Жхук1 ape. абе. д. ЛЕ Le eae. ЛЕ 
102 107 10° 10' 10° 


Work - seconds 


Figure 4.12: Work-precision diagrams for splitting methods with 256 grid points and of orders 4, 6 
and 8 applied to the NLSE. Solid lines indicate absorbing boundaries (ab) were used. 
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Figure 4.13: Work-precision diagrams for splitting methods with 512 grid points and of orders 4, 6 
and 8 applied to the NLSE. Solid lines indicate absorbing boundaries (ab) were used. 
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Figure 4.14: Work-precision diagrams for splitting methods with 1024 grid points and of orders 4, 
6 and 8 applied to the NLSE. Solid lines indicate absorbing boundaries (ab) were used. 
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error in the conservation of mass is extremely small, well below the error in the solution. One 
must, however, be careful that the step size of the method is chosen such that the error in 
the solution remains above the discretisation error. If not, the additional steps used do not 
increase the accuracy of the final solution but do cause the errors in the mass conservation 
to accumulate and increase above what they might otherwise be. 

Considering the efficiency of the splitting methods in the sense of time taken versus 
accuracy gained showed that, of the methods considered, the order four method with five 
stages is preferable for errors greater than 10-7 and that the order six method with nine 
stages is preferable for smaller errors. The methods which used additional symmetric terms 
to reduce the coefficient of the leading error term were preferable when compared with 
methods which did not. 

In the next chapter we use the techniques developed here to apply the split-step Fourier 
method to a nonlinear wave equation which describes an optical pulse propagating in a 


medium with a nonlinear amplitude response and with second order dispersion. 


5: The Split-Step Fourier Method 


for a Nonlinear Wave Equation 


5.1 Derivation of the nonlinear wave equation 


Here we derive an equation to describe the propagation of optical pulses in a nonlinear 
medium with second order dispersion. To do so we draw on the ideas introduced in chapter 


[concerning wavepackets, dispersion and nonlinear wavetrains. Further details can be found 


in [Newell and Moloney, 1992] and in |Akhmonav et al., 1992). 


We recall that in chapter[I]the dispersion relation (II) for the case of constant refractive 


index was given by 
2,2 


d(k,w)A = (e = — ) A=0. 





This was obtained by substituting a plane wave E = A exp(i(kx2—wt)) into the wave equation 


for a linear dielectric (L9) given by 


Instead of a plane wave we now insert a wavepacket of the form (L.15); F(x,t) = A(z, t) exp(i(kox— 


wot)) + c.c. where c.c. denotes the complex conjugate of the previous term. This gives 


OE 8A дА | 

m = (52 + 2iko = KA) exp(i(kox — wot)) 
ФЕ 8A .. дА | 

dE = (5 — Anz - ги) СЕИ 
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and therefore 
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which is a PDE for the slowly varying envelope of the electric field, A(x, t). 
In order to see this more directly we can explicitly consider the change in the envelope 


by calculating the derivatives of the wavepacket E(x,t) = A(x, t) exp(i(kox — wot)). 


д ) дА 
ud i(kom—wot) _ I 
Da A(x, t)e oz 
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From this we see, that for a wavepacket, the effect of applying the differential + (or 2) to 
the electric field is to multiply the envelope A(z,t) by the differential operator i (ko — ix) 
(or —i (wo + i&)). For a plane wave the effect of the same operator(s) on E(x,t) was to 
multiply the amplitude A by the algebraic coefficient ikọ (or —iwo). In both cases we can 
see that the resulting operators which we must apply to A are the necessary terms to ensure 
that the dispersion relation d(k,w) = k? — n?w?/c?=0 is satisfied. That is, for a plane wave 


k = ko, w = wo and for a wave packet 


‚0 | 
w = wo + 1. 


аат at 


We now further modify the dispersion relation to take into account the frequency de- 
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pendence of the refractive index. This is easily done by replacing the constant n with its 


frequency dependent form n(w). The dispersion relation is now 


d(k,w)A = (e E c pe 


C 


which for a wave packet gives 


А 2 А 2 
d(k,w)A = | (i = 2) Ра еза) A(z,t)=0. GD 


Ox c 


We now expand each of the terms in the dispersion relation truncating where necessary at 
order и? (recall from {Т.Т5) 2. = О(и), 25 = О(и2)): 
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To order u? the dispersion relation (B.I) is therefore given by 
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At leading order we therefore have the familiar dispersion relation kê — [n(wo)|?w2/c? = 0 


which we encountered first in section [L.1.3] 


By using the leading order dispersion relation k? — [n(wo)]?w3/c? = 0 we can rewrite the 


expression above as 
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We then use the relations 1/u = (Ik/dw)., and kz = (0?k/0u?),,, where u is known as the 
group velocity and ka as the second order dispersion coefficient, to get 
ðA | 2ko0A , 0?A ( 1 ) A 
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In the same way in which we used the leading order dispersion relation k? —[n(wo)|?w3/c? = 
0 to simplify the expansion of (B.I) we can use the dispersion relation of order и to simplify 
the order и? equation given by equation (5.2). At order и equation { Б.П) gives 
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>) A(z, t) = 0. (5.3) 
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We multiply this from the left by the formal expression — 2E MR and get 
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We сап now substitute equation (5.4) into the PDE (5.2) which allows us to eliminate the 
4 term and gives 
OA 260 OA ФА 
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a first order version of equation (6.2). Dividing through by i2ko we have 


дА 10А ik ðA _ 
euer эб С (5.5) 


An advantage of the derivation presented here is that it is easy to see how we might 
include the effect of a small perturbation ón of order и? to the refractive index. To do so 
we follow the same procedure simply replacing n(w) with n(w) + ön. When we expand the 


necessary expressions up to order u? we get 








On| ð 189m| 8? 
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(n (wo + if) + Фп)? (wo + i2) = [n (wo + үс (wo + 2) + 2wen(wo)ón + O(u?). 


We can therefore include the effect of the perturbation ón by adding an extra term 


2 2 
won(wo)dn 4 


с? 


to equation (5.2). For equation (5.5) the extra term has the form 


= tw dN 4 | iönko 


C n 


A. 


We saw in section[L.1.3]that, for a plane wave, in order to include the effect of a nonlinear 
refractive index it is necessary to replace n(w) with its intensity dependent form n(w, | A|?) = 


no(w) + n2(w)|A|?. The corresponding dispersion relation was then 


no, |A|?) Pu? 


d(k, vw, |A|?) = к? = 2 


= 0. 
c 


For a wavepacket the variation in the wave numbers and frequencies about their central 
values is of order u we therefore have |A|? = O(u?) which means that we may include the 
effect of the na] A|? term as a perturbation to the linear refractive index n(w) = no(w) of 
order и?. In other words, we can simply set ôn = n2(wo)|A|? into equations (5.2) and (5.5). 
The resulting equations therefore include the effect of the nonlinear refractive index and are 


respectively given by 





k2 
i2ko— +i — po 0 APA =0, (5.6) 
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We mention briefly one of the physical phenomena described by equations (5.6) and 
(5.7). Self-phase modulation is the change in an optical pulse due to the nonlinearity of the 
refractive index and therefore occurs for non-zero values of na. 

For an optical pulse E(x, t) = A(x, t) exp(i(kox — wot)) + c.c. with a nonlinear refractive 
index characterised by n = no + тә |А]? if the refractive index response of the medium is 
near-instantaneous then the effect of the nonlinear medium is to change the phase of the 


transmitted pulse by the amount 


—n9|A|?woL 


C 


ỌNL = 


over a propagation distance of L. If na > 0 this causes a shift to lower frequencies at 
the leading edge of the pulse while the trailing edge is shifted to higher frequencies. This 
phenomenon is known as spectral broadening ; 

In some cases (when ka < 0) this effect is counteracted by the second order (or group 
velocity) dispersion of the pulse. Balancing these two effects leads to the formation of 
optical solitons. Such pulses can travel long distances without any change in profile (as seen 
in chapter B). 

In the case when ka > 0 the effect is the opposite and both the self-phase modulation 
and the group velocity dispersion lead to dispersion of the pulse. This effect can be used as 
a simple model of the interaction of intense laser pulses with plasma where heating of the 
plasma leads to positive values of na. In general however, plasmas are highly inhomogeneous 
media and result in processes much more complicated than those described by equations 
(5.6) and (5.7); important (neglected) factors being ionisation, electron and ion density 
distributions and temporal variations (pp.548,553). Details concerning the effect 
of self-phase modulation on pulse propagation with an emphasis on the interaction with 
the group velocity dispersion and its applications to pulse propagation in optical fibres 
can be found in (pp.274-282), |Butcher and Cotter, 1990| (pp.236-245) and 
(pp.324-330, 515-517). 


5.2 Applying the split-step Fourier method 


Typically, the basic time and distance units for optical waves are of the order of femtoseconds 


(fs, or 10—15) and micrometres (um, or 1076m) respectively. Appropriate choices for the 
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parameters are: 


с = 3x 10° m/s, 
no = 1 (dimensionless), 
п € [1x10 75,1 x 10723] (dimensionless), 
ko = 7.854 x 10° mt, 


и = em. 


The effect of the (dimensionless) coefficient ka is to include second order dispersion in the 
descriptions given by equations (5.6) and (5.7). Sensible choices for the value of ka lie in the 
range ka € [0,1 x 10724]. The second order dispersion comes into play when the product 
kokz is of similar magnitude to 1/u. For simplicity, we rescale equations (5.6) and (5.7) so 
that the basic units for time and space are both unity, this allows us to drop the fs and 
um units. Setting a = 1071? and 8 = 1076, the rescaled equations for (5.6) and (5.7] are 


respectively 
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In order to model the propagation of an optical pulse suitable initial conditions are given 


by the Gaussian pulse 


A(z,t)z-0 = Ao exp (R9) = Ap exp [LÁ , (5.10) 


70 
—(t — to)? 


=. Aas (==) 2e bendi). (5.11) 
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with Ao = 3 x 1019, то = 10 and to = Зто. 


When we apply the split-step Fourier method to the problem we implicitly impose pe- 
riodic boundary conditions. In order that these do not conflict with the conditions given 
above, care must be taken to ensure that the solution remains zero in the vicinity of the 


boundaries. 


In order to apply the split-step Fourier method to equations (5.6) and (5.7) we split the 
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equations into linear and nonlinear parts. The temporal derivatives in the linear part can 
then be evaluated in Fourier space and a splitting method used to integrate with respect to 


the space variable x. 


In the case of equation (5.7) the split-step Fourier method can be applied in the same 


fashion as that described in chapter [2] For equation (5.6) we first rewrite the equation as a 








first order problem in x. We set A(x,t) = v(x, t), да = gu = u(x,t), and 
д |v 
— = 5.12 
where 
0, 1 
Е = 2 . k2 
(2) (ka 24022), —i2kopI 


We now split F as F = X + Y where the linear terms, including the time derivatives 
which are to be evaluated in Fourier space, are contained in X and the nonlinear term is 


contained in Y. That is, 





x 0 I 0 0 
= 2 : ; Y- 2 . 
(2) (kat) —i2kopI —2%0 gy 0 


We discretise along the temporal domain by setting A(z,tm) = v(zx,t,,) = V and 


94 (T, tm) = u(z,tm) = U, tm = тт +a, т = 0,...,N — 1 anda є R. We will use 
Yn to denote the discrete form of Y and x N to denote the discrete form of X in Fourier 


space. 

Because of the sparse structures of Xy and Yy, (and exp(hX м) and exp(AYy)), it is 
important that any algorithm is optimised to exploit this structure and to avoid unnecessary 
multiplications by zero. 


0 


The special structure of Yy = 
Bo 





| where Bo, is the diagonal matrix 


k2 
Ba = -2-2 g'diag(V) 
0 


means that all terms of order two and greater in the Taylor expansion of exp(hYy) are zero 
and it is possible to calculate the matrix exponential both cheaply and accurately. A direct 
calculation of exp(hX N) is more expensive. Since, for constant h, exp(hX N) is constant 


throughout the integration one possibility would be to evaluate the exponential at the start 
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of the integration and store the result for later use. An even more efficient option is to 
exploit the special structure of X N which allows use to compute the product of the matrix 
exponential with a vector exactly using trigonometric functions. 

0 

D ^I 
"y is a complex coefficient. We can therefore view 


Jem e DS I 


as a solution of the differential equation 


The matrix X N has the structure X N= | | where D is a diagonal matrix and 


with initial conditions | yo | . 
20 


In order to find an expression for | 2 | not involving a matrix exponential we solve the 
2 


system of differential equations 


у = z, (5.14) 
y" 2 = Dy+yly. (5.15) 


(g + oge”, (5.16) 
y = @ +29 + 07 Ge". (5.17) 


e 
| 


Equating (5.15) and (5.17) gives 9” = (y — 20)7 + (Doy — o?I)ğ, which after choosing 
с = y/2 becomes 
#' = (D +72/4) ӯ = —9°ў, (5.18) 
where 0 = (-D — 42 /AT)8, ie. —02 = D + y? /AI. 
In order for this definition to make sense we must ensure that 0 is positive definite. We 
have y = —2ikgoo, D = (2) (ки? - ey), и = diag (23) , n = =,...,4 — 1. 








а и 


When the parameters in these expressions are in the vicinity of those given earlier in this 


section the matrix — D — 4?/4I is positive definite for L/N many orders of magnitude larger 


106 Chapter 5. The Split-Step Fourier Method for a Nonlinear Wave Equation 


than the biggest sensible grid size for which one could hope to achieve any sort of accurate 


result. 


The transformed problem given by (5.18) can also be written in matrix notation as 


afa] fo r][à 
&[v |" [-e oll? о» 


The solution of (5.18) is given by 
g = sin (0t) Сү + cos (0t) Сэ. (5.20) 


The coefficients C1 and С» are determined by the initial conditions y(0) = %(0) and y’(0) = 
y' (0) + 0%(0) which give 


C, = 8"'(y'(0) – yy(0)/2), (5.21) 
Cy = y(0). (5.22) 


We can now calculate the product (5.13), exactly and without using any matrix-vector 
products, by using equation (5.20) and the coefficients (5.21) and (5.22) to give 


Yo 


ET e(3t), (5.23) 








(и) _ | cos(6t) 9! sin(6t) 
— | —6sin(0t) cos(bt) 








A symmetric stage of the split-step Fourier method for equation (5.6) can now be calcu- 


lated as follows: 
wt = F-lexp(hX у /2)Е exp(hYw)F-! exp(h Xy /2)FW™ (5.24) 


where 
A 0, I 


XN = 2 Й 
(2) (к? — ido) , —i2kg0I 


Ww = И" Т and u = diag (24) п = =,..., X — 1. The operators Е and F^! аге 
L 2 2 


given by 
—1 
Е Fy 0 mio Fy | 
0 Fr 0 Fx 


where Fy and Fy (= 4F%) are the FFT and IFFT operators of size N. 
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5.3 Numerical results 


We are now in a position to apply the split-step Fourier method to the first and second order 
equations, (5.8) and (6.9), for the propagation of an optical wave in a nonlinear medium. We 
make a slight modification to these equations so as to bring them in line with equations used 

у with which we will later wish to make a comparison. This 
modification consists of taking the complex conjugate of the equations and then halving the 
coefficient of the nonlinear term. Taking the complex conjugate has no effect on the physical 
implications of the equation as the electric field has the form E(z,t) = A(x, t) exp(i(kox — 
wot)) + c.c. and hence we are simply solving the PDE corresponding to the envelope of 
the second term (c.c.) of the electric field rather than the first term which we used in the 
derivation. Halving the coefficient of the nonlinear term obviously reduces the effect of this 
term which will have an impact on the physical implications of the solution. We will see in 
the section [5.3.3] the effect that this has on the final solutions. 

The modified versions of equations (5.8) and (5.9) are given respectively by 


king 
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5.3.1 Comparison of first and second order equations 


A physically sensible problem consists of considering the solution on a space interval of 100 
um and over a period of time from t = 0 to t = 300 fs. We therefore consider solutions 
defined on the scaled domain of x € [0,100] um and t € [—70, 330] fs. 

We have extended the time domain in the negative direction so as to ensure that the 
magnitude of the initial condition is close to zero near this boundary; we have A(0,—70) ~ 
10—95. We also extend the time domain slightly past the point where we wish to evaluate 
the solution so as to reduce any possible effects that the periodic boundary may have on the 
solution near the edge of the domain. 

Although all the computations involving the split-step Fourier methods in this chapter 
use the domain described here we will restrict the plots shown to the domain over which we 
are interested in the solution This is mostly x € [0,100] um, t € [0,300] fs 

We use N = 256 points along the time domain and a step size of h = 0.05 to integrate 
with respect to space. Here we take ng = 1 x 10-?4 and and kọ = 1 x 1072". All other 


parameters are the same as those given at the start of section 5.2] As we are currently 
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Figure 5.1: The solutions of the second order (upper plot) and first order (lower plot) modified wave 
equations as given by equations (5.25) and (6.26) respectively. N = 256, h = 0.05, na = 1x 10724 
and k2 = 1 x 1072". The corresponding plots for the non-modified wave equations (5.8) and (5.9) 
appear identical. 
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interested in the qualitative properties of the solutions we use the Strang splitting integrator 


which, as we saw in chapter Д performs well when high accuracy is not required. 


Figure B.J]shows the magnitudes of the solutions for the modified wave equations (6.25) 
and (5.26) over the entire domain. Because they are similar to those of the modified equations 
for this choice of parameters we omit the solutions for the original wave equations (5.8) and 
(5.9). We see immediately that while both the first and second order equations show the 
original Gaussian pulse travelling in the direction of positive t and x, the second order 
equation also shows a smaller secondary pulse with a double peak which travels in the 
direction (z,—t), that is, at the same speed as the principle solution, but in the opposite 
direction. This secondary pulse is a purely numerical phenomenon and we shall consider 


both its origins and how it may be prevented in section [5.3.2] 
We are particularly interested in the profile of the solution obtained by cutting the pulse 


in the x direction at t = 300 fs. Because we wish to use a power of two for the number of 
grid points along the time domain we have the choice of, either altering the endpoint of the 
time domain slightly (for each choice of N) so that one of the grid points lands exactly on 
t = 300 fs or, taking the grid point on the time domain closest to t = 300 fs and plotting 
the profile at that point in time. For simplicity we opt for the second of these choices. The 
profile of the non-modified version of the first order wave equation at t = 300.1961 fs is 
shown in Figure [5.2] Since (with the exception of the small, non-physical secondary peak 
mentioned above) the profiles for the other equations are similar to that of the non-modified 
first order wave equation we take the profile of the first order non-modified equation as a 
reference and subtract it from the remaining profiles and plot only the difference between 
them. This also is shown in Figure 5.2] All the profiles are taken at t = 300.1961 fs. 


We can see in FigureB.2]that for the parameters used, (ng = 1x 107-74 and ka = 1x 10-7), 
the modification to the coefficient of the nonlinear term in the wave equations has only a 
small effect. For the first order equations the difference between the modified and non- 
modified versions is almost negligible. For the second order equations the modification 
causes a change of at most about 2% in the solution. This is not the case for all choices 
of parameters, as we shall see in section [5.3.3] The small secondary peak which occurs in 
the profiles of the second order solutions at about x = 40 um is due to the secondary pulse 
travelling in the opposite direction to the principle pulse and re-entering the computational 
domain due to the periodic boundary conditions. We also see that the difference between 
the first order and second order equations is small with respect to the total amplitude of 
the solution. For both the non-modified and the modified equations the amplitude of the 
pulse for the second order equation is slightly less than that of its first order counterpart in 


the region before the peak amplitude and is slightly greater than the amplitude of the first 


110 


Chapter 5. The Split-Step Fourier Method for a Nonlinear Wave Equation 









































x 10" Profile: t = 300 fs 
3r 
2.5r 
2- 
E45. 
1+ 
0.5 - 
0 1 П L 1 | 
0 0.2 0.4 0.6 0.8 1 
х (т) x10 
x 10° Profile: t = 300 fs 
1.5 т r 
— olm - otref 
o2m - otref \ 
--- о2пт – otref | 
1- l 1 
E 
E 
B 
Ф M 
2 no^ E 
Ф 0.5; iy lh \ 7 
S ral А I | 
= I l \ ba 
ПШ ba 
voy Us 
0 j 
Y 1 
V 4 
id 
Vl 
-0.5- 1 1 
NI 
tt 
A 
=| | 1 | M 
0 0.2 0.4 0.6 0.8 1 
x (m) x 10^ 


Figure 5.2: The upper plot shows the profile of the pulse for the non-modified first order equation 
at the point on the time domain closest to t = 300 fs, (i.e. t = 300.1961 fs). Because the profiles of 
the other pulses are similar, in the lower plot we show the difference between the first order modified 
eqn. (solid line), the second order modified eqn. (dotted line) and the second order non-modified 
eqn. (dashed line) with respect to the profile for the non-modified first order eqn. (from the upper 


plot). № = 256, h = 0.05 um, no 21x 10724 and ka = 1 x 1072”. 
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order profiles afterwards. Again, we will see in section [5.3.3] that this agreement depends on 


the choice of the parameters na and ka. 


5.3.2 Characteristics and filtering for the second order wave equation 


In the plot corresponding to the second order wave equation in Figure [B.I] we saw that in 
addition to the principle pulse also visible was a secondary pulse with amplitude of slightly 
less than 296 of that of the main pulse. This secondary pulse travels towards the right of the 
domain with the same speed as the main pulse. This secondary pulse travelling in the (x, —t) 
direction occurs even though the initial conditions used correspond to a single Gaussian pulse 
travelling in the direction (x,t) This pulse is a purely numerical phenomenon and is due to 


the fact that we are solving a second order wave equation using necessarily finite step sizes. 


Characteristic curves 


'The secondary wave can be explained using the idea of characteristic curves of a PDE. For 
simplicity we consider only linear homogeneous PDEs but the same underlying idea also 


holds in the nonlinear case. If we have a first order PDE of the form 


where a and b are coefficients which possibly depend on x and t, we can parameterise т and 


t as x = z(s), t = t(s). Computing the total derivative of v with respect to s gives 


dv Ж... _ 
ds Oxds Otds | 
For de = a and u = b the above expression returns the PDE we originally considered. 
Hence, we can write us = 2 and t = т --сопі. The solutions of адо +09 = 0 can therefore 
be represented in the x,t plane by lines of gradient b/a. These are the characteristic curves 
of the PDE and are given by v(r,t) = t — Pg — const. 
If instead we consider the second order wave equation 
Әр „д% 
= - CC — =0, v=v(z,t 
Ox? ot? yo 


one finds Я = +1 and, therefore, the characteristic curves are given by «(r,t) = x + ct = 
de c 








const. 
Possible characteristic curves for both the first order and the second order case are 


illustrated in Figure [5.3] A consequence is that when the PDE for the second order system 
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i-— >т = const. 








+ cx = const. 





Figure 5.3: Possible characteristic curves in the (x,t) plane for linear homogeneous PDEs of 15% 
order (left) and 2nd order (right). 


is solved by discretising in one variable and then integrating a system of ODEs in the other 
variable, the grid-points which form the discrete domain do not, in general, match the 
characteristic curve corresponding to the initial conditions for a single pulse. The initial 
pulse is therefore effectively split into two components, one of which corresponds to a pulse 
with a characteristic curve which is approximately the same as that of the exact solution for 
the initial conditions, and a second pulse which corresponds to a characteristic curve with a 


gradient of the opposite sign to that of the original curve. 


Filtering 


The periodic boundary conditions necessary for the split-step Fourier method mean that 
if the secondary wave is allowed to propagate it will eventually leave the computational 
domain travelling in the opposite direction to that of the principle wave and will re-enter the 
domain on the opposite side. This leads to problems as the computational domain must be 
extended in the direction in which the secondary wave is propagating otherwise the solution 
for the principle wave will interact with the secondary wave introducing non-physical effects. 
We therefore wish to eliminate the secondary wave or to prevent it from re-entering the 
computational domain. To do so we can use a filter which works by multiplying the solution 
by a function which is one in the vicinity of the principle peak and is zero at the boundaries. 


Two possibilities are: 


1. A function which is zero everywhere away from the position of the principle peak such 
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Which has the form 


1 





2. A function which is zero only on the side of the domain away from the direction of 
travel of the principle peak; 





1 Xi«rz«mui 
falz) = 5 (1 + соз (=*)) dim ть (5.28) 
0 Tp < 2 < Xp. 





A possible advantage of the filter given by equation (5.28) is that it does not alter the 
solution on the section of the domain ahead of the principle peak. Although, in theory, this 
region should be close to zero anyway and hence unaffected by a filter of the form (6.27), one 
must still make a decision as to at what point on the domain the magnitude of the Gaussian 
peak is sufficiently small to be treated as zero. Using filter (5.28) removes this possibility 
but at the risk of allowing the solution to become non-zero near the boundary closest to the 
leading edge of the pulse. 

Because we are interested in the profile of the solution at t — 300 fs along the entire 
spatial domain we have used a fixed computational domain. This means that during the 
integration the pulse moves towards the boundary of the computational domain (in fact if 


we integrate to x = 100 jum it travels over the edge of the computational domain) rather 
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than remaining in the centre of the domain as it would if we used a moving domain. A filter 
of the form (5.27) is therefore poorly suited to a fixed computational domain of this size. 

For the remaining calculations we use filter (5.28) with 10% of the total domain used for 
the “dead-region” in which the solution is set to zero and 5% of the domain for the region 
in which the filter changes between zero and one. 

The fact that for an integration period of тела = 100 um part of the pulse leaves the 
computational domain means that we must extend the length of the time domain. Rather 
than adding extra points to the time domain and therefore increasing the cost of the cal- 
culation by computing the solution at points beyond those we are interested in (remember 
we are primarily interested in the solution at t = 300 fs), we instead exploit the periodic 
boundary conditions and allow the solution to reappear on the far side of the computational 
domain. 

The rationale behind this treatment is that for a sufficiently wide time domain the 
distance between the trailing edge of the original pulse and the leading edge of the pulse 
which has been produced by the periodic boundaries is great enough that the two pulses 
have no effect on each other. In order that this remain the case it is necessary that one 
is mindful of the possibility that the leading edge of the pulse may be increasing in speed 
while the trailing edge may be slowing. It is also necessary that the propagation of the pulse 
transmitted by the periodic boundaries is not altered when it appears on the far side of the 
computational domain. We therefore only apply the filter until x = 65 um. This allows us to 
prevent the secondary pulse from interfering with the solution of the principle pulse without 
altering the propagation of the principle pulse as it leaves the computational domain. 

Figure [5.4] shows the calculation for the second order wave equation repeated with the 
same parameters as used for Figure B.I] but with the filter described above. It can be seen 
from this that the filter described prevents the transmission of the unwanted secondary wave 


without affecting the propagation of the primary pulse. 


5.3.3 Effect of n and kə 


We mentioned in section [5.1] the phenomenon of self-phase modulation which can lead to 
pulse broadening for high intensity electric fields and which is due to the intensity (i.e. |A|?) 
dependence of the refractive index. This effect is included in the equations considered here 
through the presence of the coefficient na for the nonlinear term. 

A second important effect is that of second order dispersion. This effect becomes signif- 
icant when the product of the coefficient for the second order dispersion ko and the mean 


wave number kg is in the same order of magnitude as i. The dispersion resulting from 
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Figure 5.4: The solution of the modified second order wave equation (5.25) using filter (5.28). All 
parameters are the same as those used for Figure [5.1] 


non-zero values of ka is often referred to as group velocity dispersion since 


Ok ol д Ok 


Ow? дши додо 





Ко 





where u is the group velocity of the pulse. For ka > 0 the group velocity dispersion causes 
the frequency at the leading edge of the pulse to decrease while the frequency of the trailing 
edge increases. Since ok = -9u this causes the group velocity of the leading edge of 
the pulse to increase while the trailing edge slows leading to dispersion (also referred to as 
defocusing) of the pulse, and vice-versa for kg < 0. Details on this effect can be found in 
Butcher and Cotter, 1990] (pp.238-239) and (pp.277-280). 

Solutions were calculated using the modified versions of both the first and second or- 
der forms of the wave equation with the pair of coefficients (по, ко) taking the values 
(10-78, 10-27), (10-24, 10-24) and (107??, 10724). These are compared with the results from 
using (10-74, 10-2”) which are given in section 5.3.1] all other parameters remain unchanged. 

Figure [5.5] shows the magnitude of the solution for the modified second order wave 
equation (6.25) with the coefficient for the nonlinear component of the refractive index 
increased to na = 10723. All other parameters and coefficients are the same as those used 
for Figure [5.1 


We see from Figure [5.5] that as the pulse propagates it decreases in amplitude and 
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Figure 5.5: Solution of the modified second order wave equation (5.25) with na = 10778, ka = 10774. 
All other parameters are as for Figure BI] 


broadens. This broadening is most pronounced on the trailing edge of the pulse though it 
occurs also on the leading edge. The equivalent results obtained by using the modified first 
order equation (5.26) are shown in Figure [5.6] where a similar effect can be seen. In this 
case, however, the broadening appears more symmetric about the peak of the pulse. 

Since the coefficient na affects only the nonlinear term in the wave equation the effect of 
decreasing this term by a factor of a half in the modified wave equation can be inferred by 
observing the effect of changing no. 

The three-dimensional plots for the case na = 1x 10773, kə = 1x 107?" appear similar to 
those for na = 10724, ka = 1072" shown in FigureB.I]whilst those for n = 10724, kə = 1074 
are similar to those in Figures [5.5] and [5.6] but with slightly less dispersion. We therefore 
omit them. More details of the effect of na and ka on the envelope A(z,t) can be seen by 
studying the profiles of the solutions at t = 300 fs. These are shown in Figure [5.7] 

We can see from Figure[5.7]that for the second order equation, the effect of increasing n2 
to 10723 is to cause the leading edge of the pulse to steepen while for regions of larger am- 
plitude the trailing edge follows the leading edge and advances slightly. At lower amplitudes 
the profile is comparable to that for the case where ng = 10774, ko = 10727. In comparison, 
the results for the first order equation with na = 10723, ky = 1072" appear identical to those 


for по = 10-74, kə = 107?" rather than resembling the results for the second order equation 
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Figure 5.6: Solution of the modified first order wave equation (5.25) with na = 10723, kg = 10724. 
All other parameters are as for Figure [5.1] 


with ng = 10—23 and Кә = 10777. 

In order to better compare the results for the profiles in Figure [5.7] we also list the peak 
amplitude and the width for each of the profiles shown. These figures are given in Table 
[5.1] As a measure for the width we use the full-width half-maximum with respect to the 
amplitude. That is, we measure the distance between points on either side of the peak 
where the magnitude of the solution is half the value at the peak. Because the profiles 
are not always symmetric about the peak it is important that the widths in Table [Б.П are 


considered in conjunction with the profiles in Figure [5.7] 


Figure [5.7] and Table [5.T] indicate that increasing k2, the second order dispersion coeffi- 
cient causes the pulse to broaden on both the trailing and the leading edges. For the second 
order equation this broadening is more pronounced on the leading edge of the pulse while for 
the first order equation the broadening is more symmetric. In all cases, broadening of the 
pulse is accompanied by a corresponding reduction of the peak amplitude. It is interesting 
to note that although the peaks of the first and second order profiles are in good agreement 
with each other, the extent to which the pulses broaden on either side of the peaks do not 


match to the same extent. 


The difference between the results for the first and second order equations is most likely 
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Figure 5.7: Profiles of the solution for t = 300.1961 fs x € [60, 100] um for the modified first (lower 
plot) and second (upper plot) order wave equations given by (5.26) and (5.25). (note: in the lower plot 
the solid and the dotted lines coincide.) (na, k2) = (107-74, 10-27), (10—23,10—27), (10-24, 10-24), 
and (10775, 10-74). All other parameters are as for Figure 5.1] 
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na ka peak width 

(x10/0) (ит) 

order 2 3.0028 4.950 
2.9125 5.250 

1.7272 15.00 

2.1026 10.15 





order 1 | 10—24 10-?/ | 2.9993 5.000 
10-7 10-?' | 2.9931 5.000 
10-28 107-24 | 1.7214 14.85 
10-7 10-74 | 2.1030 10.10 


Table 5.1: Peak amplitudes and FWHM widths for the profiles shown in Figure [5.7] 


to have its origins in the fact that the first order equation was obtained by assuming that the 
dispersion relation of order u (6.3) was satisfied exactly and then using this expression to 
simplify the second order wave equation (5.2), thereby reducing it to a first order equation 
in space. If the order и dispersion relation (5.3) is not exact then we must expect different 
results for the first and second order forms of the wave equation. It is also possible that 
the second order wave equation is susceptible to slightly different numerical effects to its 
first order counterpart. This possibility should be further investigated by solving the second 
order problem using a different numerical method. However, as we shall see in section 5.3.4] 
other methods do not always perform well for the larger values of na and ka considered here 
making such an investigation difficult. 

As mentioned, the second order pulses tend to show the leading edge of the pulse spread- 
ing more than the trailing edge, while for the first order equations, increasing the value of ka 
tends to cause the pulse, and particularly its trailing edge, to lag slightly. For ka = 1072, 
changing n2 has only a small effect on the profiles relative to the change due to altering 
ka for the second order equation and has no noticeable effect for the first order equation. 
However, when ky = 10-74 the same change in na produces a notable change in the profiles 


for both the first and second order equations. 


5.3.4 Comparison of numerical methods 


So far all the results in this chapter have been computed using the Strang splitting method 
with a domain of x є [0,100] um, t € [-70,330] fs, an integration step size (in space) of 
x = 0.05 um, and N = 256 points for the discretisation of the time domain (i.e. a grid 
spacing of r = 1.5625 fs). We now consider the effect of using different integration methods. 


We saw in chapter Д that of the higher order splitting methods the methods of order four 
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with five stages and order six with nine stages were optimal among those considered. 

We compare the split-step Fourier methods with the Verlet method which is used with 
an FFT discretisation of the spatial domain. The resulting method is used to integrate the 
second order modified wave equation with respect to t. The Verlet method is second order 
accurate and has the advantage that it can be applied directly to the special class of ODEs 
of the form & = f(x) without first converting them to a system of first order ODEs of twice 
the dimension. 

A natural discretisation of such an ODE is given by the second order central difference 
operator 

Lnt1 — 22 + In-ı = T? f (x) (5.29) 


This three-term recurrence is the Verlet scheme in its most elementary form. To begin the 
integration with (6.29) the missing £n value can be found be introducing the quantity i = v. 


We approximate this via 
In+1l — In-1 
2T ` 


Setting n = 0 in (5.29) and (5.30) and combining the resulting equations so as to eliminate 


(5.30) 


Un = 


21] gives 
2 


E 
Ti = zo + Tto + 5 f (20). 


A more computationally useful one-step formulation of the Verlet method is given by 


first taking the approximation 
= 
Un = Un + gf (en) (5.31) 


for the value of v at the midpoint of a step. If we then proceed as we did above in order to 
eliminate х„_1 from (5.29) and (5.30) we get 


2 


Е 
Dina = 3 fan) + zd Tin: 


We now substitute (5.31) into the expression above and get 
En41 = Tn + TU, ie (5.32) 
The step is then completed by calculating v„+1 ready to use for the next step, 
T 
Un+1 = Unti + gf (n). (5.33) 


Together (5.31), (5.32) and (5.33) give a one-step implementation. When output for & = v 
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is not required (5.31) and (5.33) can be combined so that the one-step implementation of 
the Verlet method is given by 


Un 1 = Un—4 +Tf (£n), 


Intl Tn + TU ts 


Two of the numerous references which contain more details on the Verlet method are 
and [Hairer et al., 2002]. As we mentioned in the introduction, the Ver- 
let (or leap-frog) method is known to be symplectic and can also be viewed as a (Hamiltonian) 
splitting method. Here, however, we use the term “splitting methods” to refer only to those 
(operator) splitting methods already considered where the problem is split into linear and 
nonlinear parts, not to the Verlet method. 

The MatLab code for the implementation of the Verlet method used here was provided 


y [Sautter and Hochbruck, 2005 
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Figure 5.8: Profiles of the solution at t = 300 fs for the Verlet method with h = 200/512 = 0.391 шп, 
т = 0.15 fs and the Strang method with h = 0.39 um, т = 400/2560 = 0.156 fs (Strang-1) and 
h = 0.05 um, т = 400/256 ~ 1.563 fs (Strang-2). (Note: the Strang-1 and Strang-2 lines coincide.) 


We use the coefficients na = 107*4, ka = 107?" and integrate the modified second order 
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wave equation (5.25). For the Verlet method we use a spatial domain of x € [-100, 100] um. 
This ensures that the secondary pulse which occurs for this method does not affect the 
solution. The spatial domain is discretised with 512 gird points giving h = 0.390623 um. 
The Verlet method is then used to integrate from t = 0 until t = 300 fs using a step size 
of т = 0.15 fs. The profile of the solution for t = 300 fs and x € [0,100]um is shown in 
Figure [5.8] We also show the profiles as computed by the Strang splitting method; first 
using spatial and time step sizes of h = 0.39 um, т = 400/2560 > 0.156 fs, similar to those 
used for the Verlet method and then for h = 0.05 um, т = 400/256 = 1.563 fs, the step sizes 
used in the earlier calculations. The peak amplitudes and FWHM widths corresponding 
to these profiles are Apeak = 3.0526 x 1010, FWHM = 4.6875 um for the Verlet method, 
Apeak = 3.0029 x 1010, FWHM = 5.070 um for the Strang splitting with similar grid points 
to the Verlet method and Apeak = 3.0028 х 101, FWHM = 4.950 um for the Strang splitting 
with h = 0.05 um, т = 400/256 fs. We mention again briefly that the profiles for the Strang 
splitting method are evaluated at the grid point closest to t = 300 fs rather that exactly at 
this point. The difference between the point for which the profile is shown and t = 300 fs is 


at most 0.1961 fs in this case. The change in the profile due to this difference is insignificant. 


These results show two important points. Firstly, changing the grid points for the Strang 
splitting method from a fine spatial mesh and a coarse temporal mesh to a coarse spatial 
mesh and a fine temporal mesh has little effect on the solution. The profiles almost overlap, 
the peak amplitudes agree and the width of the profile with the coarse spatial mesh is only 
slightly less than that for the fine spatial mesh. The second point of interest is that both 
of the Strang splitting methods give profiles which are slightly behind the profile from the 
Verlet method. That is, the pulse for the Verlet method travels slightly faster than those 
for the Strang splitting and appears to focus slightly with the peak amplitude of the pulse 
being slightly greater than the initial amplitude. 


The most likely reason for the differences between the results of the Strang and the Verlet 
methods, and in particular the slight lag of the Strang profiles, is that the Strang methods 
use an initial condition defined along the t-axis at x = 0 while the Verlet method uses an 
initial condition defined along the x-axis for t = 0. This means that, no matter how careful 
we are in selecting the position of the initial pulse, we can not expect the solution for the 
Verlet method at x = 0 to match the initial condition for the Strang method. Two related 
factors which contribute to this are that the initial condition for the Strang method is defined 
from t = —70 fs not just from t = 0, and that each point in the initial conditions for the 
Strang method corresponds to the solution of the Verlet method after a different number of 


time steps. 


It would be possible to define the initial conditions for the Verlet method at t = —70 fs 
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rather than at t = 0 but due to the second point mentioned above this is unlikely to improve 
the match between the Verlet x = 0 solution and the Strang x = 0 initial conditions as it 
would then be necessary to take a greater number of time steps before reaching any fixed 
point in time for x = 0. Since we introduce a small error at each time step the solution 
of the Verlet method can be likened to a perturbed version of the initial condition for the 
Strang method where the perturbation increases with increasing t. 

Possible ways in which this idea could be investigated include using the Verlet method 
to integrate with respect to space (or the Strang method to integrate with respect to time) 
so that the same initial conditions can be used for both methods. Alternatively it would be 
possible to use the initial conditions from the Strang method as a boundary condition along 
x = 0 for the Verlet method with a discretisation technique such as finite differencing rather 
than the FFT method. Unfortunately, we must leave the investigation of these possibilities 


for future work. 
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Figure 5.9: Profiles at t = 300 fs showing instability of the Verlet method with n3 = 2 x 10724, 
kg = 1 x 107?" (dotted line) and n3 = 1 x 1074, kg = 2 х 107?" (solid line). All other parameters 
are the same as for Figure [5.8] 


We would also like to be able to compare the results from the Strang splitting method 
with those of the Verlet method for the larger values of na and ka considered in section [5.3.3] 
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We find, however, that with the step sizes given above for the Verlet method the method 
is unstable for values of пә and kg slightly larger than пә = 10774, & = 10777. If either 
ng or kg is increased the instability takes the form of a bulge on the computational domain 
for values of r smaller than that where the pulse is located. This bulge grows with time 
and eventually causes the solution to blow up. Figure [5.9] shows the profile of the solution 
at t = 300 fs when the Verlet method is used with (пә, о) = (2 x 10724,1 x 1072") or 
(1x 107?4,2 x 10-2"). For the same choice of the parameters na and ka the Strang splitting 
method shows a similar profile to that of the Verlet method but without the bulge and with 
the same delay observed in Figure [5.8] 


The Verlet method also fails when the number of points on the spatial domain (x € 
[—100, 100] um) is increased above N = 2332 (for 7 fixed at 0.15 fs) or when the integration 
step size is larger than about т = 0.6 fs (for N fixed at 512 points, i.e. h = 200/512 шт). 
'The upper bound on the number of points which can be used for the discretisation of the 
spatial domain with the Verlet method places a limit on the accuracy which can be obtained 
for the discretisation of this domain. We will investigate in sections [5.3.5] and [5.3.6] the 
possibility of such bounds on the spatial step size ^ and the temporal step size 7 existing 
for the splitting methods. 

If a splitting method of higher order is used instead of the Strang splitting method the 
results are qualitatively identical. Figure [5.10] shows the profiles for the solutions due to the 
order four method with five stages and the nine stage sixth order method compared with 
the profiles for the Strang and Verlet methods. 

We have seen in Figure [5.8] that altering the grid points used with the Strang splitting 
method so that the grid points match approximately those used with the Verlet method does 
not notably change the profile of the solution. We therefore use h = 0.05 wm and 400/256 fs 
for the calculations with the splitting methods in Figure [5.I0]as these give a computationally 
cheaper solution. 

Because the results for the three splitting methods are indistinguishable from one another 
when plotted in comparison to the Verlet method we take the result from the order six 
method as a reference solution and subtract it from the results for the Strang and order four 
methods. We then plot the difference between these results (also in Figure [5.10). From this 
we can see that the difference between the splitting methods is about six orders of magnitude 
smaller than the amplitude of the solution when it is greatest. We see that in this case the 
Strang method slightly underestimates the amplitude with respect to the sixth order method 
(i.e. Strang minus order six is negative) while the fourth order method gives a peak which 
is slightly larger than than of the order six solution. As one would expect the difference 


between the order four and order six methods is less than the difference between the Strang 
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Figure 5.10: The upper plot shows the profiles at t = 300 fs for the Verlet method with h = 
200/512 um, 0.15 fs (solid line) and for the Strang (dotted line), order four (dot-dash line) and order 
six (dashed line) splitting methods with h = 0.05 um, 400/256 fs. ng = 10—24, ky = 1072”. Because 
the results for the three splitting methods are indistinguishable in this plot we take the order six 
method as a reference solution and subtract it from the Strang and order four results as seen in the 
lower plot. 
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and the order six methods. 


5.3.5 Effect of changing the spatial step size 


Having seen that the Verlet method becomes unstable when the step size along the spatial 
domain is smaller than about h = 200/2382 um we now investigate to see if the same 
phenomenon is exhibited by the Strang splitting method. We used na = 10774, kg = 107?" 
and solved the same problem as that used for the Verlet method with a variety of spatial 


step sizes. A selection of representative results are shown in Figure 5.11] 
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Figure 5.11: Profiles at t ~ 300 fs for the Strang splitting method for various values of h and with 
T = 400/256 fs (solid line) and r = 400/512 fs (dotted lines). 


When the time domain is discretised with N = 256 points we find that the Strang 
splitting method becomes unstable and the solution blows up if the spatial step size is 
greater than about h = 0.4 um. In contrast to the Verlet method, however, we find no 
lower limit for the spatial step size; h = 0.001 um giving entirely stable, if computationally 
expensive, results. The Verlet method becomes unstable for a spatial step size smaller than 
h = 0.086 um. 

When step size of the Strang splitting method is decreased to 7 — 400/512 fs along the 
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time domain the instability at large step sizes appears to vanish, or rather, any instability 
only occurs when A is far larger than any usable spatial step size. We see in Figure 5.1] that 
when h = 10 um the solution with N = 512 remains stable even though the large value of 
h means that we are no longer able to accurately represent the pulse. As was the case for 
N = 256, the method remains stable as h is decreased with the solution for h = 0.001 um 
still showing no signs of instability. For clarity, and because the results for smaller h resemble 
those for h = 0.01 um we only show the profiles for h down to h = 0.01 um. 

From these observations we conclude that the Strang splitting method is not susceptible 
to the instability at small h displayed by the Verlet method. 


5.3.6 Effect of changing the temporal step size 


We saw in section [5.3.4] that with h fixed at 200/512 um the Verlet method is unstable 
when the temporal step size is greater than about 7 = 0.6 fs. This instability at large 
temporal step sizes is not as much of a concern as the instability seen for small spatial 
step sizes as in general one would like a more accurate solution and hence a smaller step 
size is necessary. However, it is possible that the instability for large step sizes can make 
calculations unnecessarily expensive by limiting the step size to below that which we would 
like to use. We therefore investigate the Strang splitting method to see if it displays similar 
behaviour. 

As we did for section [5.3.5] we take по = 10-74, ko = 107?" and solve the modified second 
order wave equation (5.25). We fix the spatial step size at h = 0.1 fs and integrate the 
equation for a range of values of т. The results are shown in Figure [5.12] 

The profiles in Figure [5.12] show that the solution remains stable over a range of values 
for т. For the case where т = 400/32 fs the solution is still stable even though the large 
temporal step size means that the number of Fourier nodes is insufficient to sample the initial 
conditions accurately and as a result, the profile no longer resembles the true solution. As 
might be expected, the results obtained using a larger value of т differ from those were т is 
smaller. 

Because the profiles are shown for the grid points closest to, rather than exactly at, 
t = 300 fs we must expect some variation in the position of the pulses shown in the profiles. 
As we are here only interested in determining whether or nor the solution becomes unstable, 
rather than giving an accurate result, this is not of concern. 

For values of r between 400/256 fs and 400/1024 fs and for 7 smaller than 400/1024 fs 
the resulting profiles are all similar in appearance so we show only two indicative examples. 
The Strang splitting method showed no signs of instability for 7 as small as 400/4096, here 
we stopped. 
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Figure 5.12: Profiles at t ~ 300 fs for the Strang splitting method with various values of r and with 
h = 0.1 um. 


The results of this section would seem to indicate that the Strang splitting method is 
less susceptible to instability from large values of the temporal step size r than the Verlet 
method. Any upper bound on the temporal step size appears to be due to the need for using 


enough points to accurately sample the solution rather than due to stability concerns. 


5.4 Conclusion 


In this chapter we have presented a derivation of a wave equation which describes the prop- 
agation of waves in a nonlinear medium. This derivation drew from the description of 
electromagnetic waves given by Maxwell’s equations as presented in the introduction. The 
resulting equations, however, can also be both derived and applied in a number of other 
contexts. In particular, early results were obtained in the study of deep water waves where 
many of the same physical phenomena seen here are also present. 

The equations developed here are presented in both first order and second order (in space) 
forms. We have seen that when the coefficients na and ka (corresponding to the nonlinear 


component of the refractive index and second order dispersion respectively) are small the 
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first and second order wave equations give comparable results. When these coefficients are 
increased the higher order effects become significant and the results of the first and second 
order equations begin to differ. 

By considering different values of na and ka we were able to illustrate the effects of self- 
phase modulation due to a nonlinear refractive index and second order or group velocity 
dispersion for a pulse travelling in a nonlinear dielectric. 

It is also possible to use a change of variables transformation to reduce the first order 
wave equation from this chapter to the nonlinear Schrödinger equation which we dealt with 
extensively in chapters [2] [3] and Д 

The results obtained using the split-step Fourier method were compared with those from 
using the Verlet method applied to the second order problem. We saw that the profiles of 
the pulses calculated using the split-step Fourier method are similar to those of the Verlet 
method but that the Verlet method shows a slight focusing effect and that the pulses of the 
split-step Fourier method lag slightly behind those of the Verlet method. It is suspected 
that this lag is due to slight differences in the initial conditions used for the two methods 
though more investigation would be needed to confirm that this is definitely the case. 

We demonstrated that the split-step Fourier method remains stable when the values of 
ng and kg are increased beyond the range for which the Verlet method is able to deliver 
solutions. This gives the split-step Fourier method an advantage over the Verlet method 
when modelling higher order effects in the wave equation. We also showed that the split- 
step Fourier method does not suffer from the instability for fine spatial grids or coarse 
temporal grids which both affect the Verlet method. 

Most of the split-step Fourier method results for the wave equations were calculated 
using the second order Strang splitting. This method was seen to be well suited to problems 
where high accuracy is not necessary. When accurate results are required we suggest the use 
of the sixth order splitting method with nine symmetric stages. This choice is based on an 
extensive study of the behaviour of a selection of symmetric splitting methods from chapter 
[4] where the splitting methods were applied to the NLSE and the results compared with an 
exact solution. 

'The performance of the methods was studied with respect to effective order behaviour, 
mass (or energy) conservation and efficiency (precision versus work). This behaviour was 
studied for splitting methods of orders one through to eight and using between one and 17 
symmetric stages. It was found that the methods of orders four and six were optimal for 
most choices of step size or desired error. All methods of higher (four or greater) order 
were capable of easily reaching an accuracy where the total error was bounded below by 


the minimum attainable spatial discretisation error which, after a quantitative study, was 
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shown to be in the order of 10711. This gave a limit for the smallest sensible step size. The 
need for using a sufficient number of grid-points along the spatial domain, (the PDE was 
discretised in space and then integrated with respect to time), was shown with the final error 
in the solution remaining bounded below, irregardless of the integration step size, when an 
insufficient number of grid-points were used. 

We saw that, although the mass conservation of the split-step Fourier methods is not 
exact, when the integration step size is chosen such that the error in the solution is not 
affected by the lower limit due to the spatial discretisation error then it is easy to obtain a 
mass conservation error not much greater than machine precision. 

The results in chapter [4] were calculated on a moving computational domain. It was 
shown that the use of such a domain does not adversely affect the accuracy of the final 
solution. 

The methods tested were all (with the exception of the first order Lie splitting method) 
symmetric methods constructed from symmetric stages and we have shown how such inte- 
grators may be derived though the use of the BCH formula and its symmetric counterpart 
to formulate a nonlinear algebraic system of equations for the coefficients necessary for the 
methods. Possible choices of the coefficients were surveyed to find those which are optimal for 
symmetric splitting methods constructed of symmetric stages. An obvious extension to this 
would be to also consider symmetric methods not necessarily constructed from symmetric 
stages of which several potential good splittings are known. The theory for the construction 
of symmetric splittings with symmetric stages (SS methods) was developed as a generalisa- 
tion of the construction of SS methods of order 2n + 2 from a composition of methods of 
order 2n, n € N. Such methods have easily determined exact coefficients but are both more 
expensive and less accurate then the SS methods of fewer stages which were later developed. 

The results in chapter ]show that there is a benefit to be gained by using extra stages in 
a method so as to provide free variables which can be used to ‘tune’ the method by reducing 
the coefficient of the leading error term. It was found that this additional accuracy is in 
the range of one and a half to half an order or magnitude in the final error for the problem 
considered. The results in chapter Yalso illustrated the importance of considering the effect 
of the boundary conditions on non-periodic solutions as it was seen that using an absorbing 
potential has a negative effect on the mass conservation properties of the method. 

In chapter [2]we developed the general idea of a splitting method based on ideas presented 
in the introduction. We described how such a method can be applied to a nonlinear PDE 
after discretising the PDE so as to obtain a system of ODEs. We then showed how geometric 
properties of the splitting method can be used to preserve geometric properties present in the 


original problem, in particular, conservation of mass. The example of the NLSE was used to 
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illustrate conservation, dispersion and instability properties of the analytic problem. We then 
looked at how these same properties are affected by the discretisation of the PDE in the cases 
of both finite differencing and Fourier discretisations. Numerical results were presented for 
the case of periodic plane wave solutions of the NLSE. These results confirmed the theory 
that it is possible to reproduce the analytic instability of the 1-D NLSE and illustrated 
how the preservation of the mass conservation property of the NLSE can lead to the well 
known recurrence phenomenon. Theoretical step size restrictions for preventing numerical 
instability were derived and verified numerically for both spatial discretisation methods. It 
was also found that even when using an optimised version of Gaussian elimination to exploit 
the special structure of the problem, the finite difference method was less efficient than the 
Fourier method. 

In chapter B]the NLSE results from chapter [2] were extended (for the Fourier method) to 
the case of non-periodic solutions — the “fundamental solutions” of nonlinear problems. It 
was shown that the split-step Fourier method can be used both cheaply and easily to give 
numerical soliton solutions to the NLSE and that such solutions are able to approximately 
preserve many of the important physical properties of the system; in particular, mass con- 
servation and the description of a collision between solitons travelling in the same medium. 
The split-step Fourier method can therefore be used along with the NLSE (or the nonlinear 
wave equation which it can be derived from) to provide a good description of many of the 


important physical properties of waves propagating in nonlinear media. 


5.4.1 Future work 


In addition to illustrating some of the properties of the nonlinear wave equation and the 
results that one may expect when solving it using a split-step Fourier method we have also 
raised a number of possibilities for future work. 

We saw that (for the parameters at which a comparison could be made) the split-step 
Fourier methods produced results which did not entirely match those from the Verlet method 
which was used to give a comparison. To see why this is the case one should first ensure 
that the initial conditions for the two methods are equivalent. This could mean using the 
Verlet method to integrate with respect to space (as the split-step Fourier method does in 
this case) so that the same initial conditions could be used for both methods. (Alternatively, 
the split-step Fourier method could be used to integrate with respect to time and use the 
initial conditions of the Verlet method.) 

If the Verlet method was used with a finite difference discretisation of the spatial domain 
rather than the FF'T discretisation used here, then it would be possible to use the initial 


conditions from the split-step Fourier method as one of the boundary conditions for the 


132 Chapter 5. The Split-Step Fourier Method for a Nonlinear Wave Equation 


Verlet method when integrating with respect to time. This could provide another interesting 
comparison for the solution of the split-step Fourier method. 

The solution of the nonlinear wave equations should also be further investigated for the 
case where na and ka are large. To do this it would be good to have a numerical method 
other than the split-step Fourier method which could be used as a comparison. Instability 
rules out the Verlet method from such a role. Alternatively, a better understanding of the 
solutions of the split-step Fourier method for when na and ka are large could come from 
further analytic study of the nonlinear wave equations and the physical phenomena which 
they represent. 

Related to this is the study of the nonlinear wave equations when the amplitude of the 
pulse is altered. For large values of Ao the Verlet method fails and the split-step Fourier 
method gives solutions where the pulse shows lots of dispersion on the edge where z is largest. 
Because the coefficients involved in such a calculation span many orders of magnitude it 
would be useful to study which effects are in danger of being lost or altered due to finite (in 
this case 32 bit) precision. This could be investigated numerically by using a code which 
employs higher (for example, 64 bit) precision. 

It would also be interesting to extend the investigation to include other splitting methods. 
All the splitting methods considered here are symmetric methods composed of symmetric 
stages. The types of methods considered could be extended to include symmetric methods 
not composed of symmetric stages and also methods which are not symmetric at all, i.e, of 


the form 
k 


Ф(А) = П exp(a;AA) exp(b;AB) 
i=1 

where the coefficients a; and b; need satisfy no conditions other than those necessary for 
the method to be of a particular order. Such methods should be considered with respect 
to their accuracy, efficiency and treatment of various physical/geometric properties of the 
differential equations of interest. 

Beyond nonlinear oscillatory problems and wave-type equations there lies a wide field 
of problems where splitting methods can be investigated and tested. These range from 
molecular dynamics to quantum and quantum statistical mechanics to particle accelerator 


physics and celestial mechanics. 
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